# LDS project: replication data
# Clean up
rm(list = ls())
# Install Packages
library(foreign)
library(haven)
library(stargazer)
library(xtable)
library(readxl)
library(lmtest)
library(lattice)
library(dplyr)
library(multiwayvcov)
library(sandwich)
library(ggplot2)  
library(CBPS)
library(factoextra)
library(RItools)
library(psych)
library(ltm)
library(ggpubr)
library(ggmap)
library(choroplethr)
library(choroplethrMaps)
library(USAboundaries) # For a set of US map data points. 
library(USAboundariesData) # Datasets related to USAboundaries
library(ggplot2)
library(tidyverse)
library(tigris)
library(effsize)

# set working directory
setwd()
# Read Data
# lds survey data frame
# current version "lds_replication.rds"
lds_replication_data <- readRDS("lds_replication.rds")
lds_replication <- data.frame(lds_replication_data)
# geographical information for figure 1 usa map
full_lds_usa_data <- readRDS("full_lds_usa.rds")
full_lds_usa <- data.frame(full_lds_usa_data)
# geographical information for figure 1 world map
full_lds_world_data <- readRDS("full_lds_world.rds")
full_lds_world <- data.frame(full_lds_world_data)
# missionary assignment: country
lds_country_count_data <- readRDS("lds_country_count.rds")
lds_country_count <- data.frame(lds_country_count_data)
# missionary assignment: us states
lds_usa_state_count_data <- readRDS("lds_usa_state_count.rds")
lds_usa_state_count <- data.frame(lds_usa_state_count_data)
# foreign-born population by state
acs_foreign_born_state_2010_data <- readRDS("acs_foreign_born_state_2010.rds")
acs_foreign_born_state_2010 <- data.frame(acs_foreign_born_state_2010_data)
# latino population by state
acs_latino_state_2010_data <- readRDS("acs_latino_state_2010.rds")
acs_latino_state_2010 <- data.frame(acs_latino_state_2010_data)
# state-level measurement of immigration attitudes
cces_state_immi_est_data <- readRDS("/cces_state_immi_est.rds")
cces_state_immi_est <- data.frame(cces_state_immi_est_data)
# cces 2014-2014 panel
cces_panel_2012_2014_data <- readRDS("cces_panel_2012_2014.rds")
cces_panel_2012_2014 <- data.frame(cces_panel_2012_2014_data)
# wvs survey
wvs_data <- readRDS("wvs.rds")
wvs <- data.frame(wvs_data)
# international migrant
intl_migrant_data <- readRDS("intl_migrant.rds")
intl_migrant <- data.frame(intl_migrant_data)
# country, US states, and US county descriptive statistics
country_state_county_desc_data <- readRDS("country_state_county_desc.rds")
country_state_county_desc <- data.frame(country_state_county_desc_data)
# geographical information for creating the map
# us states
geo_states <- readRDS("geo_states.rds")
# world
geo_world <- readRDS("geo_world.rds")
# register google api
register_google(key = "AIzaSyCl26aKcJyjB4fAnWGU0Tx33wZBushVAX8")

########replication code for "How Social Context Affects Immigration Attitudes"########
## Table 1: Demographics of LDS Missionaries (Pre-Mission Survey)
# gender
# proportion of female missionaries
mean_female <- mean(lds_replication$female, na.rm = TRUE)
# standard deviation of female
sd_female <- sd(lds_replication$female, na.rm = TRUE)
# mean age at the time of signing up
mean_age <- mean(lds_replication$age, na.rm = TRUE)
# standard deviation of age at the time of signing up
sd_age <- sd(lds_replication$age, na.rm = TRUE)
# proportion of white missionaries
mean_white <- mean(lds_replication$white, na.rm = TRUE)
# standard deviation white missionaries
sd_white <- sd(lds_replication$white, na.rm = TRUE)
# mean family income
mean_family_income <- mean(lds_replication$family_income, na.rm = TRUE)
# standard deviation family income
sd_family_income <- sd(lds_replication$family_income, na.rm = TRUE)
# proportion of missionaries had prior foreign language experiences
mean_prior_other_lang <- mean(lds_replication$prior_other_lang, na.rm = TRUE)
# standard deviation prior foreign language experiences
sd_prior_other_lang <- sd(lds_replication$prior_other_lang, na.rm = TRUE)
# proportion of missionaries had prior health issues
mean_prior_health_issue <- mean(lds_replication$prior_health_issue, na.rm = TRUE)
# standard deviation prior health issues
sd_prior_health_issue <- sd(lds_replication$prior_health_issue, na.rm = TRUE)
# mean party id
mean_pid_republican <- mean(lds_replication$pid_republican, na.rm = TRUE)
# standard deviation party id
sd_pid_republican <- sd(lds_replication$pid_republican, na.rm = TRUE)
# mean ideology
mean_ideo_conservative <- mean(lds_replication$ideo_conservative, na.rm = TRUE)
# standard deviation ideology
sd_ideo_conservative <- sd(lds_replication$ideo_conservative, na.rm = TRUE)
# combine data
# vector collecting means
vec_desc_means <- c(mean_female, mean_age, mean_white, mean_family_income, mean_prior_other_lang, mean_prior_health_issue, mean_pid_republican, mean_ideo_conservative)
vec_desc_sds <- c(sd_female, sd_age, sd_white, sd_family_income, sd_prior_other_lang, sd_prior_health_issue, sd_pid_republican, sd_ideo_conservative)
# create table
mat_desc_demog <- cbind(vec_desc_means, vec_desc_sds)
dat_desc_demog <- as.data.frame(mat_desc_demog)
# column names
colnames(dat_desc_demog) <- c("Mean", "Standard Deviation")
# row names
row.names(dat_desc_demog) <- c("Female", "Age", "White", "Family Income", "Prior Foreign Language", "Prior Health Issue", "Republican", "Cosnervative")
# make Table 1: Demographics of LDS Missionaries (Pre-Mission Survey)
tab_desc_demog <- xtable(dat_desc_demog, digits = 2, caption = "Descriptive Statistics of the lds_replication Missionaries
(Pre-Mission Survey)")
print(tab_desc_demog)

## statistics on page 11
# months between post-mission survey completion and returning from mission
mean(lds_replication$diff_month_complete, na.rm = TRUE)
# mean approx. 6 months
median(lds_replication$diff_month_complete, na.rm = TRUE)
# median = 3 months
sd(lds_replication$diff_month_complete, na.rm = TRUE)
# sd approx.y 7 months

## Figure 1: Mission Service Locations
# USA map
plot_mission_location_usa <- ggplot() + 
  geom_sf(data = full_lds_usa, aes(geometry = geometry, fill = scaled_count), size = 0.04, color = "grey") +
  geom_sf(data = geo_states, aes(fill = NA), size = 0.2, colour = "#333333") + 
  scale_fill_brewer(palette = "YlOrRd",
                    labels = c("1", "2", "3", "4", "5", "5-10", "10+", "None")) + 
  theme_void(base_size = 16) + 
  labs(fill = "") + 
  theme(legend.key.size = unit(0.3, 'cm'),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 8),
        legend.position = c(0.93, 0.5)) + 
  coord_sf()
# world map
plot_mission_location_world <- ggplot() + 
  geom_map(data = full_lds_world, map = geo_world,
           aes(long, lat, map_id = country_name, fill = scaled_count), 
           size = 0.2, color = "black") + 
  scale_fill_brewer(palette = "YlOrRd",
                    labels = c("1-5", "5-10", "10-15", "15-20", "20-25", "25-30", "30+", "None")) + 
  theme_void() + 
  labs(fill = "") + 
  theme(legend.key.size = unit(0.3, 'cm'),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 8),
        legend.position = c(0.9, 1.1)) + 
  coord_fixed()
# combine two maps
bind_usa_world_map <- ggarrange(plot_mission_location_usa, 
                                plot_mission_location_world,
                                widths = c(1, 1),
                                ncol = 1, 
                                nrow = 2,
                                common.legend = FALSE, 
                                align = "v",
                                legend = "right")
# annotate the combined plot
## Figure 1: Mission Service Locations
plot_bind_usa_world_map <- annotate_figure(bind_usa_world_map,
                                           top = text_grob("LDS Missionaries Around the World", 
                                                           color = "Black", 
                                                           face = "bold", 
                                                           size = 20))
print(plot_bind_usa_world_map)

## statistics on page 13
print(table(lds_replication$usa, exclude = NULL)[1] / (table(lds_replication$usa, exclude = NULL)[1] + table(lds_replication$usa, exclude = NULL)[2]))
# approx. 65% served abroad

## create Table 2A: Effects of Larger Latino/Immigrant Communities: Policy Items, US Sample
# usa sample
lds_replication_usa <- subset(lds_replication, lds_replication$usa == 1)
dim(lds_replication_usa)
# abroad sample
lds_replication_abroad <- subset(lds_replication, lds_replication$usa == 0)
dim(lds_replication_abroad)
# vector containing demographic variable name
var_demog <- c("female", "age", "white", "family_income", "prior_health_issue", "prior_other_lang", "pid_republican", "ideo_conservative")
# variable names for regression tables
var_name_demog <- c("Female", "Age", "White", "Family Income", "Prior Health Issue", "Prior Foreign Language", "Republican", "Conservative")
# outcome variables
var_diff_immi_outcome <- c("diff_immi_policy", "diff_immi_stereotype")
# latino, policy items
f_diff_county_latino_policy <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_loc_county_latino", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_latino_policy_mod <- lm(formula = f_diff_county_latino_policy, data = lds_replication_usa)
summary(diff_county_latino_policy_mod)
# foreign-born, policy items
f_diff_county_foreign_born_policy <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_loc_county_foreign_born", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_foreign_born_policy_mod <- lm(formula = f_diff_county_foreign_born_policy, data = lds_replication_usa)
summary(diff_county_foreign_born_policy_mod)
# foreign-born from latin american countries, policy items
f_diff_county_fbla_policy <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_loc_county_fbla", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_fbla_policy_mod <- lm(formula = f_diff_county_fbla_policy, data = lds_replication_usa)
summary(diff_county_fbla_policy_mod)
## Table 2A: Effects of Larger Latino/Immigrant Communities: Policy Items, US Sample
# variable names
var_name_county_diff <- c("Diff. in % Latino", "Diff. in % Foreign Born", "Diff in % Foreign Born from Latin America")
var_name_county_diff_policy_outcome <- c("Policy Items")
# regression table
stargazer(diff_county_latino_policy_mod, 
          diff_county_foreign_born_policy_mod, 
          diff_county_fbla_policy_mod,
          omit = c(var_demog,
                   "mission_year"),
          title = "Effects of Larger Latino/Immigrant Communities: Stereotype Items, US Sample",
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          font.size = "small",
          covariate.labels = var_name_county_diff, 
          dep.var.labels = var_name_county_diff_policy_outcome,
          add.lines = list(c("Demographics?", "Yes", "Yes", "Yes"), c("Mission Year FE?", "Yes", "Yes", "Yes")),
          label = "out_county_diff_policy")

## create Table 2B: Effects of Larger Latino/Immigrant Communities: Stereotype Items, US Sample
# latino, stereotype items
f_diff_county_latino_stereotype <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_loc_county_latino", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_latino_stereotype_mod <- lm(formula = f_diff_county_latino_stereotype, data = lds_replication_usa)
summary(diff_county_latino_stereotype_mod)
# foreign-born, stereotype items
f_diff_county_foreign_born_stereotype <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_loc_county_foreign_born", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_foreign_born_stereotype_mod <- lm(formula = f_diff_county_foreign_born_stereotype, data = lds_replication_usa)
summary(diff_county_foreign_born_stereotype_mod)
# foreign-born from latin american countrie, stereotype items
f_diff_county_fbla_stereotype <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_loc_county_fbla", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_fbla_stereotype_mod <- lm(formula = f_diff_county_fbla_stereotype, data = lds_replication_usa)
summary(diff_county_fbla_stereotype_mod)
# create Table 2B: Effects of Larger Latino/Immigrant Communities: stereotype Items, US Sample
# variable names
var_name_county_diff_stereotype_outcome <- c("Stereotype Items")
# regression table
stargazer(diff_county_latino_stereotype_mod, 
          diff_county_foreign_born_stereotype_mod, 
          diff_county_fbla_stereotype_mod,
          omit = c(var_demog,
                   "mission_year"),
          title = "Effects of Larger Latino/Immigrant Communities: Stereotype Items, US Sample",
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          font.size = "small",
          covariate.labels = var_name_county_diff, 
          dep.var.labels = var_name_county_diff_stereotype_outcome,
          add.lines = list(c("Demographics?", "Yes", "Yes", "Yes"), c("Mission Year FE?", "Yes", "Yes", "Yes")),
          label = "out_county_diff_stereotype_mod")

## Figure 2: Foreign-Born-Latin American Community and Immigration Attitudes
# scatter plot for policy items
plot_diff_fbla_policy <- ggplot() +
  geom_point(data = lds_replication_usa, 
             aes(x = diff_loc_county_fbla, 
                 y = diff_immi_policy),
             size = 1) +
  geom_smooth(data = lds_replication_usa, 
              aes(x = diff_loc_county_fbla, 
                  y = diff_immi_policy),
              method = "lm", 
              formula = y ~ x, se = FALSE) +
  ylab(NULL) +
  xlab(NULL) +
  ylim(c(-1, 1)) +
  xlim(c(-50, 50)) +
  ggtitle("Policy", subtitle = NULL) +
  scale_colour_brewer(palette = "Set1", direction = -1) +
  theme_light(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5))
# scatter plot for stereotype items
plot_diff_fbla_stereotype <- ggplot() +
  geom_point(data = lds_replication_usa, 
             aes(x = diff_loc_county_fbla, 
                 y = diff_immi_stereotype),
             size = 1) +
  geom_smooth(data = lds_replication_usa, 
              aes(x = diff_loc_county_fbla, 
                  y = diff_immi_stereotype),
              method = "lm", 
              formula = y ~ x, se = FALSE) +
  ylab(NULL) +
  xlab(NULL) +
  ylim(c(-1, 1)) +
  xlim(c(-50, 50)) +
  ggtitle("Stereotype", subtitle = NULL) +
  scale_colour_brewer(palette = "Set1", direction = -1) +
  theme_light(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5))
# combine two plots
bind_diff_fbla <- ggarrange(plot_diff_fbla_policy,
                            plot_diff_fbla_stereotype,
                            widths = c(0.3, 0.3),
                            heights = c(0.5, 0.5), 
                            ncol = 2, 
                            nrow = 1,
                            align = "hv",
                            legend = "right",
                            common.legend = TRUE)
# annotate
## Figure 2: Foreign-Born-Latin-American Community and Immigration Attitudes
plot_bind_diff_fbla <- annotate_figure(bind_diff_fbla,
                                       top = text_grob("Foreign-Born-Latin-American Community and \n Immigration Attitudes", 
                                                       color = "Black", 
                                                       face = "bold", 
                                                       size = 14),
                                       bottom = text_grob("% Foreign-Born From Latin America in Mission Service Locations - % Foreign-Born From Latin America in Pre-Mission Most-Identified Places", 
                                                          color = "Black", 
                                                          face = "bold", 
                                                          size = 10),
                                       left = text_grob("Change in Immigration Attitudes", color = "Black", face = "bold", size = 10, rot = 90))
print(plot_bind_diff_fbla)

## Table 3: Effects of Immigration Attitudes of the Local People: US Sample
# policy items
# with foreign-born as covariate only
f_diff_local_attitude_policy_1 <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_local_attitude", "mission_loc_county_foreign_born"), collapse = " + "), sep = " ~ "))
diff_local_attitude_policy_mod_1 <- lm(formula = f_diff_local_attitude_policy_1, data = lds_replication_usa)
summary(diff_local_attitude_policy_mod_1)
# with covariates
f_diff_local_attitude_policy_2 <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_local_attitude", "mission_loc_county_foreign_born", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_local_attitude_policy_mod_2 <- lm(formula = f_diff_local_attitude_policy_2, data = lds_replication_usa)
summary(diff_local_attitude_policy_mod_2)
# stereotype items
# with foreign-born as covariate only
f_diff_local_attitude_stereotype_1 <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_local_attitude", "mission_loc_county_foreign_born"), collapse = " + "), sep = " ~ "))
diff_local_attitude_stereotype_mod_1 <- lm(formula = f_diff_local_attitude_stereotype_1, data = lds_replication_usa)
summary(diff_local_attitude_stereotype_mod_1)
# with covariates
f_diff_local_attitude_stereotype_2 <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_local_attitude", "mission_loc_county_foreign_born", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_local_attitude_stereotype_mod_2 <- lm(formula = f_diff_local_attitude_stereotype_2, data = lds_replication_usa)
summary(diff_local_attitude_stereotype_mod_2)
## Table 3: Effects of Immigration Attitudes of the Local People: US Sample
# variable names
var_name_local_attitude_outcome <- c("Policy", "Stereotype")
var_name_local_attitude <- c("Attitudes of the Locals", "% Foreign Born", var_name_demog)
# regression table
stargazer(diff_local_attitude_policy_mod_1,
          diff_local_attitude_policy_mod_2,
          diff_local_attitude_stereotype_mod_1, 
          diff_local_attitude_stereotype_mod_2,
          omit = c("mission_year"),
          title = "Effects of Immigration Attitudes of the Local People: US Sample",
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_local_attitude, 
          dep.var.labels = var_name_local_attitude_outcome,
          font.size = "small",
          add.lines = list(c("Mission Year FE?", "No", "Yes", "No", "Yes")),
          label = "diff_local_attitude")

## create Figure 3: Effects of being Assigned to Speak a Foreign Language, USA Sample
# run the regression model
# policy items
# with covariates
f_non_eng_mission_usa_policy <- formula(paste(var_diff_immi_outcome[1], paste(c("non_eng_mission", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_usa_policy_mod <- lm(formula = f_non_eng_mission_usa_policy, data = lds_replication_usa)
summary(non_eng_mission_usa_policy_mod)
# stereotype items
# with covariates
f_non_eng_mission_usa_stereotype <- formula(paste(var_diff_immi_outcome[2], paste(c("non_eng_mission", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_usa_stereotype_mod <- lm(formula = f_non_eng_mission_usa_stereotype, data = lds_replication_usa)
summary(non_eng_mission_usa_stereotype_mod)
# extract coefficients to create plot
# policy items 
# point estimate
coef_non_eng_mission_usa_policy <- summary(non_eng_mission_usa_policy_mod)$coefficients[2, 1]
# standard error
se_non_eng_mission_usa_policy<- summary(non_eng_mission_usa_policy_mod)$coefficients[2, 2]
# combine
qoi_non_eng_mission_usa_policy <- as.vector(c(coef_non_eng_mission_usa_policy, se_non_eng_mission_usa_policy))
# confidence intervals
ci_non_eng_mission_usa_policy <- as.vector(confint(non_eng_mission_usa_policy_mod)[2, ])
# add confidence intervals
qoi_non_eng_mission_usa_policy <- c(qoi_non_eng_mission_usa_policy, 
                                    ci_non_eng_mission_usa_policy)
qoi_non_eng_mission_usa_policy <- as.data.frame(qoi_non_eng_mission_usa_policy)
# readjust
qoi_non_eng_mission_usa_policy <- t(qoi_non_eng_mission_usa_policy)
qoi_non_eng_mission_usa_policy <- as.data.frame(qoi_non_eng_mission_usa_policy)
# column names
colnames(qoi_non_eng_mission_usa_policy) <- c("point_est", 
                                              "se", 
                                              "non_eng_mission_usa_lower_95", 
                                              "non_eng_mission_usa_upper_95")
# clean row names
row.names(qoi_non_eng_mission_usa_policy) <- c()
# extract coefficients to create plot
# stereotype items 
# point estimate
coef_non_eng_mission_usa_stereotype <- summary(non_eng_mission_usa_stereotype_mod)$coefficients[2, 1]
# standard error
se_non_eng_mission_usa_stereotype<- summary(non_eng_mission_usa_stereotype_mod)$coefficients[2, 2]
# combine
qoi_non_eng_mission_usa_stereotype <- as.vector(c(coef_non_eng_mission_usa_stereotype, se_non_eng_mission_usa_stereotype))
# confidence intervals
ci_non_eng_mission_usa_stereotype <- as.vector(confint(non_eng_mission_usa_stereotype_mod)[2, ])
# add confidence intervals
qoi_non_eng_mission_usa_stereotype <- c(qoi_non_eng_mission_usa_stereotype, 
                                        ci_non_eng_mission_usa_stereotype)
qoi_non_eng_mission_usa_stereotype <- as.data.frame(qoi_non_eng_mission_usa_stereotype)
# readjust
qoi_non_eng_mission_usa_stereotype <- t(qoi_non_eng_mission_usa_stereotype)
qoi_non_eng_mission_usa_stereotype <- as.data.frame(qoi_non_eng_mission_usa_stereotype)
# column names
colnames(qoi_non_eng_mission_usa_stereotype) <- c("point_est", 
                                                  "se", 
                                                  "non_eng_mission_usa_lower_95", 
                                                  "non_eng_mission_usa_upper_95")
# clean row names
row.names(qoi_non_eng_mission_usa_stereotype) <- c()
# combine policy and stereotype results
out_non_eng_mission_usa <- rbind(qoi_non_eng_mission_usa_policy, qoi_non_eng_mission_usa_stereotype)
out_non_eng_mission_usa <- as.data.frame(out_non_eng_mission_usa)
# add policy domain indicator
out_non_eng_mission_usa$Domain <- rep(NA)
out_non_eng_mission_usa$Domain <- c("Policy", "Stereotype")
## Figure 3: Effects of being Assigned to Speak a Foreign Language, USA Sample
plot_non_eng_mission_usa <- ggplot() +
  geom_pointrange(data = out_non_eng_mission_usa, 
                  aes(x = Domain, 
                      y = point_est, 
                      ymin = non_eng_mission_usa_lower_95, 
                      ymax = non_eng_mission_usa_upper_95,
                      shape = Domain),
                  position = position_dodge(width = 0.4),
                  size = 0.8) +
  annotate("text",
           x = c(2.1, 0.9),
           y = c(-0.055, -0.1),
           label = as.character(rev(round(out_non_eng_mission_usa$point_est, 3)))) +  ylab("Change in Immigration Attitudes") +
  xlab("Issue Domain") +
  geom_hline(yintercept = 0, linetype="dashed") +
  ylim(c(-0.2, 0.2)) +
  ggtitle("Effects of being Assigned to Speak a Foreign Language, USA Sample") +
  scale_colour_brewer(palette = "Set1", direction=-1) +
  theme_light(base_size = 14) +
  coord_flip() +
  theme(text = element_text(size = 16), 
        plot.title = element_text(hjust = 0.5))
print(plot_non_eng_mission_usa)

## statistics on page 26 Cohen��s D for the effect of the assignment to speak a foreign language
# policy items
# speak foreign language
non_eng_mission_usa_diff_immi_policy <- lds_replication_usa$diff_immi_policy[lds_replication_usa$non_eng_mission == 1]
# speak english
eng_mission_usa_diff_immi_policy <- lds_replication_usa$diff_immi_policy[lds_replication_usa$non_eng_mission == 0]
# drop missing observations
non_eng_mission_usa_diff_immi_policy <- non_eng_mission_usa_diff_immi_policy[complete.cases(non_eng_mission_usa_diff_immi_policy)]
eng_mission_usa_diff_immi_policy <- eng_mission_usa_diff_immi_policy[complete.cases(eng_mission_usa_diff_immi_policy)]
# compute Cohen��s D for policy items
cohen_d_non_eng_mission_usa_policy <-  cohen.d(non_eng_mission_usa_diff_immi_policy, 
                                               eng_mission_usa_diff_immi_policy)
print(cohen_d_non_eng_mission_usa_policy$estimate)
# Cohen��s D for policy items = 0.57
# stereotype items
# speak foreign language
non_eng_mission_usa_diff_immi_stereotype <- lds_replication_usa$diff_immi_stereotype[lds_replication_usa$non_eng_mission == 1]
# speak english
eng_mission_usa_diff_immi_stereotype <- lds_replication_usa$diff_immi_stereotype[lds_replication_usa$non_eng_mission == 0]
# drop missing observations
non_eng_mission_usa_diff_immi_stereotype <- non_eng_mission_usa_diff_immi_stereotype[complete.cases(non_eng_mission_usa_diff_immi_stereotype)]
eng_mission_usa_diff_immi_stereotype <- eng_mission_usa_diff_immi_stereotype[complete.cases(eng_mission_usa_diff_immi_stereotype)]
# compute Cohen��s D for stereotype items
cohen_d_non_eng_mission_usa_stereotype <-  cohen.d(non_eng_mission_usa_diff_immi_stereotype, 
                                               eng_mission_usa_diff_immi_stereotype)
print(cohen_d_non_eng_mission_usa_stereotype$estimate)
# Cohen��s D for stereotype items = 0.20

## create Table 4: Effects of Larger Immigrant Communities
# policy items
# no covariate
f_intl_migrant_policy_1 <- formula(paste(var_diff_immi_outcome[1], c("loc_country_intl_migrant"), sep = " ~ "))
intl_migrant_policy_mod_1 <- lm(formula = f_intl_migrant_policy_1, data = lds_replication)
summary(intl_migrant_policy_mod_1)
# with covariates
f_intl_migrant_policy_2 <- formula(paste(var_diff_immi_outcome[1], paste(c("loc_country_intl_migrant", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
intl_migrant_policy_mod_2 <- lm(formula = f_intl_migrant_policy_2, data = lds_replication)
summary(intl_migrant_policy_mod_2)
# stereotype items
# no covariate
f_intl_migrant_stereotype_1 <- formula(paste(var_diff_immi_outcome[2], c("loc_country_intl_migrant"), sep = " ~ "))
intl_migrant_stereotype_mod_1 <- lm(formula = f_intl_migrant_stereotype_1, data = lds_replication)
summary(intl_migrant_stereotype_mod_1)
# with covariates
f_intl_migrant_stereotype_2 <- formula(paste(var_diff_immi_outcome[2], paste(c("loc_country_intl_migrant", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
intl_migrant_stereotype_mod_2 <- lm(formula = f_intl_migrant_stereotype_2, data = lds_replication)
summary(intl_migrant_stereotype_mod_2)
## Table 4: Effects of Larger Immigrant Communities (% International Migrant)
# variable names
var_name_intl_migrant_outcome <- c("Policy", "Stereotype")
var_name_intl_migrant <- c("% International Migrant")
# regression table
stargazer(intl_migrant_policy_mod_1,
          intl_migrant_policy_mod_2,
          intl_migrant_stereotype_mod_1,
          intl_migrant_stereotype_mod_2,
          title = "Effects of Larger Immigrant Communities (% International Migrant)",
          omit = c(var_demog, "mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_intl_migrant, 
          dep.var.labels = var_name_intl_migrant_outcome,
          font.size = "small",
          add.lines = list(c("Demographics", "No", "Yes", "No", "Yes"), 
                           c("Year FE?", "No", "Yes", "No", "Yes")),
          label = "diff_imm_intl_migrant")

## Appendix
## make Table A1: Summary of International Missionary Locations (Countries)
# select columns
out_lds_country_count <- lds_country_count[c("country_name", "count")]
# column names
colnames(out_lds_country_count) <- c("Country", "Count")
## Table A1: Summary of International Missionary Locations (Countries)
# table output
tab_lds_country_count <- xtable(out_lds_country_count, digits = 0, caption = "Summary of International Missionary Locations (Countries)")
print(tab_lds_country_count)

## make Table A2: Summary of Missionary Locations (US States)
# data for output
out_lds_usa_state_count <- as.data.frame(lds_usa_state_count)
# column names
colnames(out_lds_usa_state_count) <- c("State", "Count")
## Table A2: Summary of Missionary Locations (US States)
# table output
tab_lds_usa_state_count <- xtable(out_lds_usa_state_count, digits = 0, caption = "Summary of Missionary Locations (US States)")
print(tab_lds_usa_state_count)

## make Table A3: Balance Tests: Served in USA vs. Abroad, Demographics
# specify variables of interest
var_x_bal_usa_abroad_full <- c(var_demog, "usa")
# select variables
dat_x_bal_usa_abroad_full <- lds_replication[var_x_bal_usa_abroad_full]
# drop missing observations
dat_x_bal_usa_abroad_full <- dat_x_bal_usa_abroad_full[complete.cases(dat_x_bal_usa_abroad_full), ]
# specify formula
f_x_bal_usa_abroad_full <- formula(paste("usa", paste(var_demog, collapse = " + "), sep = " ~ "))
# xbalance function: usa vs. abroad, demographics
out_x_bal_usa_abroad_full <- xBalance(f_x_bal_usa_abroad_full,
                                      data = dat_x_bal_usa_abroad_full,
                                      na.rm = TRUE,
                                      # specify statistics to present
                                      report = c("adj.means", 
                                                 "adj.mean.diffs",
                                                 "std.diffs",
                                                 "z.scores", 
                                                 "p.values"))
out_x_bal_usa_abroad_full <- as.data.frame(out_x_bal_usa_abroad_full)
# column names
colnames(out_x_bal_usa_abroad_full) <- c("Abroad", "USA", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_usa_abroad_full) <- var_name_demog
## Table A3: Balance Tests: Served in USA vs. Abroad, Demographics
## NOTE: I swapped the location of prior foreign language and prior health issue when presenting the balance table ##
# create latex table
tab_x_bal_usa_abroad_full <- xtable(out_x_bal_usa_abroad_full, digits = 2, caption = "Balance Tests: Served in USA vs. Abroad, Demographics")
print(tab_x_bal_usa_abroad_full)
# check n for usa and abroad
table(dat_x_bal_usa_abroad_full$usa, exclude = NULL)
# N_USA = 902; N_Abroad = 477
# stratified by gender
# specify formula
# drop female as a covariate here
f_x_bal_usa_abroad_full_str_gender <- formula(paste("usa", paste(var_demog[! (var_demog %in% c("female"))], collapse = " + "), sep = " ~ ")) 
# xbalance function: usa vs. abroad, demographics, stratified by gender
out_x_bal_usa_abroad_full_str_gender <- xBalance(f_x_bal_usa_abroad_full_str_gender,
                               strata = factor(dat_x_bal_usa_abroad_full$female),
                               data = dat_x_bal_usa_abroad_full,
                               na.rm = TRUE,
                               # specify statistics to present
                               report = c("adj.means", 
                                          "adj.mean.diffs", 
                                          "std.diffs",
                                          "z.scores", 
                                          "p.values"))
out_x_bal_usa_abroad_full_str_gender <- as.data.frame(out_x_bal_usa_abroad_full_str_gender)
# column names
colnames(out_x_bal_usa_abroad_full_str_gender) <- c("Abroad", "USA", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_usa_abroad_full_str_gender) <- var_name_demog[! (var_name_demog %in% c("Female"))]
## Table A3: Balance Tests: Served in USA vs. Abroad, Demographics (Stratified by Gender)
# create latex table
tab_x_bal_usa_abroad_full_str_gender <- xtable(out_x_bal_usa_abroad_full_str_gender, digits = 2, caption = "Balance Tests: Served in USA vs. Abroad, Demographics (Stratified by Gender)")
print(tab_x_bal_usa_abroad_full_str_gender)

## make Table A4: Balance Test: Served in USA vs. Abroad, Pre-Mission Immigration Attitudes
# specify variables of interest
# pre-mission immigration attitude
var_immi_attitude_pre <- c("immi_policy_pre",
                           "immi_stereotype_pre",
                           "immi_english_pre",
                           "immi_tuition_pre",
                           "immi_welfare_pre",
                           "immi_employer_pre",
                           "immi_legal_pre",
                           "immi_rights_pre",
                           "immi_violence_pre",
                           "immi_lazy_pre",
                           "immi_honest_pre")
# post-mission immigration attitude
var_immi_attitude_post <- c("immi_policy_post",
                            "immi_stereotype_post",
                            "immi_english_post",
                            "immi_tuition_post",
                            "immi_welfare_post",
                            "immi_employer_post",
                            "immi_legal_post",
                            "immi_rights_post",
                            "immi_violence_post",
                            "immi_lazy_post",
                            "immi_honest_post")
# add usa-abroad indicator
var_x_bal_immi_attitude_usa_abroad <- c("usa", var_immi_attitude_pre)
# select variables
dat_x_bal_immi_usa_abroad_full <- lds_replication[var_x_bal_immi_attitude_usa_abroad]
# drop missing observations
dat_x_bal_immi_usa_abroad_full <- dat_x_bal_immi_usa_abroad_full[complete.cases(dat_x_bal_immi_usa_abroad_full), ]
# specify formula
f_x_bal_immi_usa_abroad_full <- formula(paste("usa", paste(var_immi_attitude_pre, collapse = " + "), sep = " ~ "))
# xbalance function: usa vs. abroad, pre-mission immigration attitudes
out_x_bal_immi_usa_abroad_full <- xBalance(f_x_bal_immi_usa_abroad_full,
                                           data = dat_x_bal_immi_usa_abroad_full,
                                           na.rm = TRUE,
                                           # specify statistics to present
                                           report = c("adj.means", 
                                                      "adj.mean.diffs",
                                                      "std.diffs",
                                                      "z.scores", 
                                                      "p.values"))
out_x_bal_immi_usa_abroad_full <- as.data.frame(out_x_bal_immi_usa_abroad_full)
# column names
colnames(out_x_bal_immi_usa_abroad_full) <- c("Abroad", "USA", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_immi_usa_abroad_full) <- c("Policy", "Stereotype", "English", "Tuition", "Welfare", "Employer", "Legal Status", "Rights", "Violence", "Lazy", "Honest")
## Table A4: Balance Test: Served in USA vs. Abroad, Pre-Mission Immigration Attitudes
# create latex table
tab_x_bal_immi_usa_abroad_full <- xtable(out_x_bal_immi_usa_abroad_full, digits = 2, caption = "Balance Test: Served in USA vs. Abroad, Pre-Mission Immigration Attitudes")
print(tab_x_bal_immi_usa_abroad_full)
# check n for usa and abroad
table(dat_x_bal_immi_usa_abroad_full$usa, exclude = NULL)
# N_USA = 984; N_Abroad = 539

## make Table A5: Balance Test: Pre-Mission Attitudes and Assigned to Speak Spanish, USA Sample
# add spanish mission indicator
var_x_bal_immi_attitude_esp_mission_usa <- c("esp_mission", var_immi_attitude_pre)
# select variables
dat_x_bal_immi_esp_mission_usa <- lds_replication_usa[var_x_bal_immi_attitude_esp_mission_usa]
# drop missing observations
dat_x_bal_immi_esp_mission_usa <- dat_x_bal_immi_esp_mission_usa[complete.cases(dat_x_bal_immi_esp_mission_usa), ]
# specify formula
f_x_bal_immi_esp_mission_usa <- formula(paste("esp_mission", paste(var_immi_attitude_pre, collapse = " + "), sep = " ~ "))
# xbalance function: assigned to speak spanish, pre-mission immigration attitudes 
out_x_bal_immi_esp_mission_usa <- xBalance(f_x_bal_immi_esp_mission_usa,
                                           data = dat_x_bal_immi_esp_mission_usa,
                                           na.rm = TRUE,
                                           # specify statistics to present
                                           report = c("adj.means", 
                                                      "adj.mean.diffs",
                                                      "std.diffs",
                                                      "z.scores", 
                                                      "p.values"))
out_x_bal_immi_esp_mission_usa <- as.data.frame(out_x_bal_immi_esp_mission_usa)
# column names
colnames(out_x_bal_immi_esp_mission_usa) <- c("Other Languages", "Speak Spanish", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_immi_esp_mission_usa) <- c("Policy", "Stereotype", "English", "Tuition", "Welfare", "Employer", "Legal Status", "Rights", "Violence", "Lazy", "Honest")
## Table A5: Balance Test: Pre-Mission Attitudes and Assigned to Speak Spanish, USA Sample
# create latex table
tab_x_bal_immi_esp_mission_usa <- xtable(out_x_bal_immi_esp_mission_usa, digits = 2, caption = "Balance Test: Pre-Mission Attitudes and Assigned to Speak Spanish, USA Sample")
print(tab_x_bal_immi_esp_mission_usa)
# check n for usa and abroad
table(dat_x_bal_immi_esp_mission_usa$esp_mission, exclude = NULL)
# N_Spanish = 113; N_Other_Languages = 360

## make Table A6: Balance Test: Demographics and Assigned to Speak Spanish, USA Sample
###############################NOTE#############################################
## we used those for whom we have full data (pre- and post-mission            ## 
## immigration questions) to generate the results for table a6.               ##      
## If we add those with missing data (pre-mission attitudes only respondents),##
## the results are nearly identical (tiny substantive changes) but that prior ##
##language experience reaches statistical significance.                       ##
################################NOTE############################################
# specify variables of interest
var_x_bal_demog_esp_mission_usa <- c(var_x_bal_immi_attitude_esp_mission_usa, var_demog, var_immi_attitude_post)
# select variables
dat_x_bal_demog_esp_mission_usa <- lds_replication_usa[var_x_bal_demog_esp_mission_usa]
# drop missing observations
dat_x_bal_demog_esp_mission_usa <- dat_x_bal_demog_esp_mission_usa[complete.cases(dat_x_bal_demog_esp_mission_usa), ]
# specify formula
f_x_bal_demog_esp_mission_usa <- formula(paste("esp_mission", paste(var_demog, collapse = " + "), sep = " ~ "))
# xbalance function: assigned to speak spanish and demographics, usa sample
out_x_bal_demog_esp_mission_usa <- xBalance(f_x_bal_demog_esp_mission_usa,
                                            data = dat_x_bal_demog_esp_mission_usa,
                                            na.rm = TRUE,
                                            # specify statistics to present
                                            report = c("adj.means", 
                                                       "adj.mean.diffs",
                                                       "std.diffs",
                                                       "z.scores", 
                                                       "p.values"))
out_x_bal_demog_esp_mission_usa <- as.data.frame(out_x_bal_demog_esp_mission_usa)
# column names
colnames(out_x_bal_demog_esp_mission_usa) <- c("Other Languages", "Speak Spanish", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_demog_esp_mission_usa) <- var_name_demog
## Table A6: Balance Test: Demographics and Assigned to Speak Spanish, USA Sample
# create latex table
tab_x_bal_demog_esp_mission_usa <- xtable(out_x_bal_demog_esp_mission_usa, digits = 2, caption = "Balance Test: Pre-Mission Attitudes and Assigned to Speak Spanish, USA Sample (Drop Missing Post-Attitudes)")
print(tab_x_bal_demog_esp_mission_usa)
# check n for usa and abroad
table(dat_x_bal_demog_esp_mission_usa$esp_mission, exclude = NULL)
# N_Spanish = 96; N_Other_Languages = 310

## make Table A7: Balance Test: Pre-Mission Attitudes and Assigned to Speak Foreign Language, USA Sample
# specify variables of interest
var_x_bal_immi_attitude_non_eng_mission_usa <- c("non_eng_mission", var_immi_attitude_pre)
# select variables
dat_x_bal_immi_attitude_non_eng_mission_usa <- lds_replication_usa[var_x_bal_immi_attitude_non_eng_mission_usa]
# drop missing observations
dat_x_bal_immi_attitude_non_eng_mission_usa <- dat_x_bal_immi_attitude_non_eng_mission_usa[complete.cases(dat_x_bal_immi_attitude_non_eng_mission_usa), ]
# specify formula
f_x_bal_immi_attitude_non_eng_mission_usa <- formula(paste("non_eng_mission", paste(var_immi_attitude_pre, collapse = " + "), sep = " ~ "))
# xbalance function: assigned to speak spanish and demographics, usa sample
out_x_bal_immi_attitude_non_eng_mission_usa <- xBalance(f_x_bal_immi_attitude_non_eng_mission_usa,
                                                        data = dat_x_bal_immi_attitude_non_eng_mission_usa,
                                                        na.rm = TRUE,
                                                        # specify statistics to present
                                                        report = c("adj.means", 
                                                                   "adj.mean.diffs",
                                                                   "std.diffs",
                                                                   "z.scores", 
                                                                   "p.values"))
out_x_bal_immi_attitude_non_eng_mission_usa <- as.data.frame(out_x_bal_immi_attitude_non_eng_mission_usa)
# column names
colnames(out_x_bal_immi_attitude_non_eng_mission_usa) <- c("Speak English", "Speak Other Language", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_immi_attitude_non_eng_mission_usa) <- c("Policy", "Stereotype", "English", "Tuition", "Welfare", "Employer", "Legal Status", "Rights", "Violence", "Lazy", "Honest")
## Table A7: Balance Test: Pre-Mission Attitudes and Assigned to Speak Foreign Language, USA Sample
# create latex table
tab_x_bal_immi_attitude_non_eng_mission_usa <- xtable(out_x_bal_immi_attitude_non_eng_mission_usa, digits = 2, caption = "Balance Test: Pre-Mission Attitudes and Assigned to Speak Foreign Language, USA Sample")
print(tab_x_bal_immi_attitude_non_eng_mission_usa)
# check n for usa and abroad
table(dat_x_bal_immi_attitude_non_eng_mission_usa$non_eng_mission, exclude = NULL)
# N_Non_English = 130; N_English = 343

## make Table A8: Demographics and Survey Completion, Full Sample
# specify variables of interest
var_x_bal_demog_completion_full <- c("complete", var_demog)
# select variables
dat_x_bal_demog_completion_full <- lds_replication[var_x_bal_demog_completion_full]
# drop missing observations
dat_x_bal_demog_completion_full <- dat_x_bal_demog_completion_full[complete.cases(dat_x_bal_demog_completion_full), ]
# specify formula
f_x_bal_demog_completion_full <- formula(paste("complete", paste(var_demog, collapse = " + "), sep = " ~ "))
# xbalance function: assigned to speak spanish and demographics, usa sample
out_x_bal_demog_completion_full <- xBalance(f_x_bal_demog_completion_full,
                                            data = dat_x_bal_demog_completion_full,
                                            na.rm = TRUE,
                                            # specify statistics to present
                                            report = c("adj.means", 
                                                       "adj.mean.diffs",
                                                       "std.diffs",
                                                       "z.scores", 
                                                       "p.values"))
out_x_bal_demog_completion_full <- as.data.frame(out_x_bal_demog_completion_full)
# column names
colnames(out_x_bal_demog_completion_full) <- c("Incomplete", "Complete", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_demog_completion_full) <- var_name_demog
## Table A8: Demographics and Survey Completion, Full Sample
# create latex table
tab_x_bal_demog_completion_full <- xtable(out_x_bal_demog_completion_full, digits = 2, caption = "Balance Test: Demographics and Survey Completion, Full Sample")
print(tab_x_bal_demog_completion_full)
# check n for usa and abroad
table(dat_x_bal_demog_completion_full$complete, exclude = NULL)
# N_Complete = 1332; N_Incomplete = 173

## make Table A9: Balance Test: Pre-Mission Attitudes and Survey Completion, Full Sample
# specify variables of interest
var_x_bal_immi_attitude_completion_full <- c("complete", var_immi_attitude_pre)
# select variables
dat_x_bal_immi_attitude_completion_full <- lds_replication[var_x_bal_immi_attitude_completion_full]
# drop missing observations
dat_x_bal_immi_attitude_completion_full <- dat_x_bal_immi_attitude_completion_full[complete.cases(dat_x_bal_immi_attitude_completion_full), ]
# specify formula
f_x_bal_immi_attitude_completion_full <- formula(paste("complete", paste(var_immi_attitude_pre, collapse = " + "), sep = " ~ "))
# xbalance function: assigned to speak spanish and demographics, usa sample
out_x_bal_immi_attitude_completion_full <- xBalance(f_x_bal_immi_attitude_completion_full,
                                            data = dat_x_bal_immi_attitude_completion_full,
                                            na.rm = TRUE,
                                            # specify statistics to present
                                            report = c("adj.means", 
                                                       "adj.mean.diffs",
                                                       "std.diffs",
                                                       "z.scores", 
                                                       "p.values"))
out_x_bal_immi_attitude_completion_full <- as.data.frame(out_x_bal_immi_attitude_completion_full)
# column names
colnames(out_x_bal_immi_attitude_completion_full) <- c("Incomplete", "Complete", "Adjusted-Diff", "Std-Diff", "z", "p-value")
##NOTE: swapped the location of employer and legal status##
# row names
row.names(out_x_bal_immi_attitude_completion_full) <- c("Policy", "Stereotype", "English", "Tuition", "Welfare", "Employer", "Legal Status", "Rights", "Violence", "Lazy", "Honest")
## Table A9: Balance Test: Pre-Mission Attitudes and Survey Completion, Full Sample
# create latex table
tab_x_bal_immi_attitude_completion_full <- xtable(out_x_bal_immi_attitude_completion_full, digits = 2, caption = "Balance Test: Demographics and Survey Completion, Full Sample")
print(tab_x_bal_immi_attitude_completion_full)
# check n for usa and abroad
table(dat_x_bal_immi_attitude_completion_full$complete, exclude = NULL)
# N_Complete = 1466; N_Incomplete = 198

## make Table A10: Balance Test: Most Identified Places
# size of immigrant community in most-idnetified places
var_immi_demog_identified_place <- c("identified_county_latino", "identified_county_foreign_born", "identified_county_fbla")
# specify variables of interest
var_x_bal_usa_abroad_immi_demog_identified_place_full <- c("usa", var_immi_demog_identified_place)
# select variables
dat_x_bal_usa_abroad_immi_demog_identified_place_full <- lds_replication[var_x_bal_usa_abroad_immi_demog_identified_place_full]
# drop missing observations
dat_x_bal_usa_abroad_immi_demog_identified_place_full <- dat_x_bal_usa_abroad_immi_demog_identified_place_full[complete.cases(dat_x_bal_usa_abroad_immi_demog_identified_place_full), ]
# specify formula
f_x_bal_usa_abroad_immi_demog_identified_place_full <- formula(paste("usa", paste(var_immi_demog_identified_place, collapse = " + "), sep = " ~ "))
# xbalance function: usa vs. abroad and identified places
out_x_bal_usa_abroad_immi_demog_identified_place_full <- xBalance(f_x_bal_usa_abroad_immi_demog_identified_place_full,
                                                                  data = dat_x_bal_usa_abroad_immi_demog_identified_place_full,
                                                                  na.rm = TRUE,
                                                                  # specify statistics to present
                                                                  report = c("adj.means", 
                                                                             "adj.mean.diffs",
                                                                             "std.diffs",
                                                                             "z.scores", 
                                                                             "p.values"))
out_x_bal_usa_abroad_immi_demog_identified_place_full <- as.data.frame(out_x_bal_usa_abroad_immi_demog_identified_place_full)
# column names
colnames(out_x_bal_usa_abroad_immi_demog_identified_place_full) <- c("Abroad", "USA", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_usa_abroad_immi_demog_identified_place_full) <- c("% Latino in Most-Identified Places", "% Foreign-Born in Most-Identified Places", "% Foreign-Born from Latin America in Most-Identified Places")
## Table A10(1): Balance Test: USA vs. Abroad and Most Identified Places
# create latex table
tab_x_bal_usa_abroad_immi_demog_identified_place_full <- xtable(out_x_bal_usa_abroad_immi_demog_identified_place_full, digits = 2, caption = "Balance Test: Most Identified Places, Full Sample")
print(tab_x_bal_usa_abroad_immi_demog_identified_place_full)
# check n for usa and abroad
table(dat_x_bal_usa_abroad_immi_demog_identified_place_full$usa, exclude = NULL)
# N_Abroad = 921; N_USA = 503
## make Table A10(2): Balance Test: USA vs. Abroad and Most Identified Places, Full Sample, (Stratified by Gender)
# stratified by gender
# specify variables of interest
var_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender <- c("usa", "female", var_immi_demog_identified_place)
# select variables
dat_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender <- lds_replication[var_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender]
# drop missing observations
dat_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender <- dat_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender[complete.cases(dat_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender), ]
# specify formula
f_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender <- formula(paste("usa", paste(var_immi_demog_identified_place, collapse = " + "), sep = " ~ ")) 
# xbalance function: usa vs. abroad, demographics, stratified by gender
out_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender <- xBalance(f_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender,
                                                                             strata = factor(dat_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender$female),
                                                                             data = dat_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender,
                                                                             na.rm = TRUE,
                                                                             # specify statistics to present
                                                                             report = c("adj.means", 
                                                                             "adj.mean.diffs", 
                                                                             "std.diffs",
                                                                             "z.scores", 
                                                                             "p.values"))
out_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender <- as.data.frame(out_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender)
# column names
colnames(out_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender) <- c("Abroad", "USA", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender) <- c("% Latino in Most-Identified Places", "% Foreign-Born in Most-Identified Places", "% Foreign-Born from Latin America in Most-Identified Places")
## Table A10(2): Balance Test: USA vs. Abroad and Most Identified Places (Stratified by Gender)
# create latex table
tab_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender <- xtable(out_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender, digits = 2, caption = "Balance Test: USA vs. Abroad and Most Identified Places, Full Sample, (Stratified by Gender)")
print(tab_x_bal_usa_abroad_immi_demog_identified_place_full_str_gender)
## make Table A10(3): Balance Test: Assigned to Speak Spanish and Most Identified Places, USA Sample
# specify variables of interest
var_x_bal_esp_mission_immi_demog_identified_place_usa <- c("esp_mission", var_immi_demog_identified_place)
# select variables
dat_x_bal_esp_mission_immi_demog_identified_place_usa <- lds_replication_usa[var_x_bal_esp_mission_immi_demog_identified_place_usa]
# drop missing observations
dat_x_bal_esp_mission_immi_demog_identified_place_usa <- dat_x_bal_esp_mission_immi_demog_identified_place_usa[complete.cases(dat_x_bal_esp_mission_immi_demog_identified_place_usa), ]
# specify formula
f_x_bal_esp_mission_immi_demog_identified_place_usa <- formula(paste("esp_mission", paste(var_immi_demog_identified_place, collapse = " + "), sep = " ~ "))
# xbalance function: assigned to speak spanish and identified places
out_x_bal_esp_mission_immi_demog_identified_place_usa <- xBalance(f_x_bal_esp_mission_immi_demog_identified_place_usa,
                                                                  data = dat_x_bal_esp_mission_immi_demog_identified_place_usa,
                                                                  na.rm = TRUE,
                                                                  # specify statistics to present
                                                                  report = c("adj.means", 
                                                                             "adj.mean.diffs",
                                                                             "std.diffs",
                                                                             "z.scores", 
                                                                             "p.values"))
out_x_bal_esp_mission_immi_demog_identified_place_usa <- as.data.frame(out_x_bal_esp_mission_immi_demog_identified_place_usa)
# column names
colnames(out_x_bal_esp_mission_immi_demog_identified_place_usa) <- c("Other Language", "Speak Spanish", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_esp_mission_immi_demog_identified_place_usa) <- c("% Latino in Most-Identified Places", "% Foreign-Born in Most-Identified Places", "% Foreign-Born from Latin America in Most-Identified Places")
## Table A10(3): Balance Test: Assigned to Speak Spanish and Most Identified Places, USA Sample
# create latex table
tab_x_bal_esp_mission_immi_demog_identified_place_usa <- xtable(out_x_bal_esp_mission_immi_demog_identified_place_usa, digits = 2, caption = "Assigned to Speak Spanish and Most Identified Places, USA Sample")
print(tab_x_bal_esp_mission_immi_demog_identified_place_usa)
# check n for assigned to speak spanish/other languages
table(dat_x_bal_esp_mission_immi_demog_identified_place_usa$esp_mission, exclude = NULL)
# N_Spanish = 107; N_Other_Language = 338
## make Table A10(4): Balance Test: Survey Completion and Most Identified Places, Full Sample
# specify variables of interest
var_x_bal_completion_immi_demog_identified_place_full <- c("complete", var_immi_demog_identified_place)
# select variables
dat_x_bal_completion_immi_demog_identified_place_full <- lds_replication[var_x_bal_completion_immi_demog_identified_place_full]
# drop missing observations
dat_x_bal_completion_immi_demog_identified_place_full <- dat_x_bal_completion_immi_demog_identified_place_full[complete.cases(dat_x_bal_completion_immi_demog_identified_place_full), ]
# specify formula
f_x_bal_completion_immi_demog_identified_place_full <- formula(paste("complete", paste(var_immi_demog_identified_place, collapse = " + "), sep = " ~ "))
# xbalance function: usa vs. abroad and identified places
out_x_bal_completion_immi_demog_identified_place_full <- xBalance(f_x_bal_completion_immi_demog_identified_place_full,
                                                                  data = dat_x_bal_completion_immi_demog_identified_place_full,
                                                                  na.rm = TRUE,
                                                                  # specify statistics to present
                                                                  report = c("adj.means", 
                                                                             "adj.mean.diffs",
                                                                             "std.diffs",
                                                                             "z.scores", 
                                                                             "p.values"))
out_x_bal_completion_immi_demog_identified_place_full <- as.data.frame(out_x_bal_completion_immi_demog_identified_place_full)
# column names
colnames(out_x_bal_completion_immi_demog_identified_place_full) <- c("Incomplete", "Complete", "Adjusted-Diff", "Std-Diff", "z", "p-value")
# row names
row.names(out_x_bal_completion_immi_demog_identified_place_full) <- c("% Latino in Most-Identified Places", "% Foreign-Born in Most-Identified Places", "% Foreign-Born from Latin America in Most-Identified Places")
## Table A10(4): Balance Test: Survey Completion and Most Identified Places, Full Sample
# create latex table
tab_x_bal_completion_immi_demog_identified_place_full <- xtable(out_x_bal_completion_immi_demog_identified_place_full, digits = 2, caption = "Balance Test: Survey Completion and Most Identified Places, Full Sample")
print(tab_x_bal_completion_immi_demog_identified_place_full)
# check n for complete and incomplete
table(dat_x_bal_completion_immi_demog_identified_place_full$complete, exclude = NULL)
# N_Complete = 1317; N_Incomplete = 245

## make Table A11: Balance Check: Latino, Foreign-Born, and Foreign-Born from Latin America Population, and Covariates
# outcome variables
var_bal_mission_loc_immi <- c("mission_loc_county_latino", "mission_loc_county_foreign_born", "mission_loc_county_fbla")
# latino
# specify formula
f_bal_mission_loc_latino <- formula(paste(var_bal_mission_loc_immi[1], paste(c("immi_policy_pre", "immi_stereotype_pre", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
bal_mission_loc_latino_mod <- lm(formula = f_bal_mission_loc_latino, data = lds_replication_usa)
summary(bal_mission_loc_latino_mod)
# foreign-born
# specify formula
f_bal_mission_loc_foreign_born <- formula(paste(var_bal_mission_loc_immi[2], paste(c("immi_policy_pre", "immi_stereotype_pre", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
bal_mission_loc_foreign_born_mod <- lm(formula = f_bal_mission_loc_foreign_born, data = lds_replication_usa)
summary(bal_mission_loc_foreign_born_mod)
# foreign-born from latin american countries
# specify formula
f_bal_mission_loc_fbla <- formula(paste(var_bal_mission_loc_immi[3], paste(c("immi_policy_pre", "immi_stereotype_pre", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
bal_mission_loc_fbla_mod <- lm(formula = f_bal_mission_loc_fbla, data = lds_replication_usa)
summary(bal_mission_loc_fbla_mod)
## Table A11: Balance Check: Latino, Foreign-Born, and Foreign-Born from Latin America Population, and Covariates
# variable names
var_name_bal_mission_loc_immi <- c("Pre-Mission Policy Attitudes", "Pre-Mission Stereotype Attitudes", var_name_demog)
var_name_bal_mission_loc_immi_outcome <- c("% Latino", "% Foreign Born", "% Foreign Born from Latin America")
# regression table
stargazer(bal_mission_loc_latino_mod, 
          bal_mission_loc_foreign_born_mod, 
          bal_mission_loc_fbla_mod,
          omit = c("mission_year"),
          title = "Balance Check: Latino, Foreign-Born, and Foreign-Born from Latin America Population, and Covariates",
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          font.size = "small",
          covariate.labels = var_name_bal_mission_loc_immi, 
          dep.var.labels = var_name_bal_mission_loc_immi_outcome,
          add.lines = list(c("Mission Year FE?", "Yes", "Yes", "Yes")),
          label = "out_bal_mission_loc_mod")

## make Table A12: Balance Check: Mission Service Locations and Most-Identified Pre-Mission Places, USA Sample
# latino
# specify formula
f_immi_demog_identified_place_latino <- formula(paste(var_bal_mission_loc_immi[1], paste(c(var_immi_demog_identified_place, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
bal_immi_demog_identified_place_latino_mod <- lm(formula = f_immi_demog_identified_place_latino, data = lds_replication_usa)
summary(bal_immi_demog_identified_place_latino_mod)
# foreign-born
# specify formula
f_immi_demog_identified_place_foreign_born <- formula(paste(var_bal_mission_loc_immi[2], paste(c(var_immi_demog_identified_place, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
bal_immi_demog_identified_place_foreign_born_mod <- lm(formula = f_immi_demog_identified_place_foreign_born, data = lds_replication_usa)
summary(bal_immi_demog_identified_place_foreign_born_mod)
# foreign-born from Latin American couintries
# specify formula
f_immi_demog_identified_place_fbla <- formula(paste(var_bal_mission_loc_immi[3], paste(c(var_immi_demog_identified_place, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
bal_immi_demog_identified_place_fbla_mod <- lm(formula = f_immi_demog_identified_place_fbla, data = lds_replication_usa)
summary(bal_immi_demog_identified_place_fbla_mod)
## Table A12: Balance Check: Mission Service Locations and Most-Identified Pre-Mission Places, USA Sample
# variable names
var_name_bal_immi_demog_identified_place <- c("% Latino in Most-Identified Places", "% Foreign-Born in Most-Identified Places", "% Foreign-Born from Latin America in Most-Identified Places")
var_name_bal_immi_demog_identified_place_outcome <- as.vector(var_name_bal_mission_loc_immi_outcome)
# regression table
stargazer(bal_immi_demog_identified_place_latino_mod, 
          bal_immi_demog_identified_place_foreign_born_mod, 
          bal_immi_demog_identified_place_fbla_mod,
          omit = c("mission_year"),
          title = "Balance Check: Latino, Foreign-Born, and Foreign-Born from Latin America Population, and Covariates",
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          font.size = "small",
          covariate.labels = var_name_bal_immi_demog_identified_place, 
          dep.var.labels = var_name_bal_immi_demog_identified_place_outcome,
          add.lines = list(c("Mission Year FE?", "Yes", "Yes", "Yes")),
          label = "out_bal_immi_demog_identified_placec_mod")

## make Table A13: Distribution of States with which Missionaries Most Identified Pre-Mission
# get the frequency of most-identified states
tab_loc_identified_state <- as.data.frame(table(lds_replication$loc_identified_state, exclude = NULL))
# us states only
# get the vector of us state abbreviations
state_abb <- as.vector(state.abb)
tab_loc_identified_state <- subset(tab_loc_identified_state, tab_loc_identified_state$Var1 %in% state_abb)
colnames(tab_loc_identified_state) <- c("State Abbreviations", "Frequency")
# re-order by frequency
tab_loc_identified_state_ordered <- tab_loc_identified_state[order(-tab_loc_identified_state$Frequency), ]
row.names(tab_loc_identified_state_ordered) <- c()
# Table A13: Distribution of States with which Missionaries Most Identified Pre-Mission
# latex outputs
out_tab_loc_identified_state_ordered <- xtable(tab_loc_identified_state_ordered, digits = 0, caption = "Distribution of States with which Missionaries Most Identified Pre-Mission")
print(out_tab_loc_identified_state_ordered)

## make Table A14: Table A14: Distribution of Cities in Utah with which Missionaries Most Identified Pre-Mission
# Utah identified subset
lds_identified_ut <- subset(lds_replication, lds_replication$loc_identified_state == "UT")
tab_lds_identified_ut_city <- as.data.frame(table(lds_identified_ut$loc_identified_city, exclude = NULL))
# column names
colnames(tab_lds_identified_ut_city) <- c("City", "Count")
# re-order by frequency
tab_lds_identified_ut_city_ordered <- tab_lds_identified_ut_city[order(-tab_lds_identified_ut_city$Count), ]
row.names(tab_lds_identified_ut_city_ordered) <- c()
# drop count of NAs
tab_lds_identified_ut_city_ordered <- subset(tab_lds_identified_ut_city_ordered, is.na(tab_lds_identified_ut_city_ordered$City) == FALSE)
## Table A14: Table A14: Distribution of Cities in Utah with which Missionaries Most Identified Pre-Mission
# latex outputs
out_tab_lds_identified_ut_city_ordered <- xtable(tab_lds_identified_ut_city_ordered, digits = 0, caption = "Distribution of Cities in Utah with which Missionaries Most Identified Pre-Mission")

## make Table A15: Demographics and Most Identified Places
# outcome variables
var_bal_identified_loc <- c("provo", "orem", "salt_lake_city", "utah")
# Provo
# specify formula
f_bal_identified_provo <- formula(paste(var_bal_identified_loc[1], paste(c(var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
# logit model
bal_identified_provo_mod <- glm(formula = f_bal_identified_provo, data = lds_replication, family=binomial(link = "logit"))
summary(bal_identified_provo_mod)
# Orem
# specify formula
f_bal_identified_orem <- formula(paste(var_bal_identified_loc[2], paste(c(var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
# logit model
bal_identified_orem_mod <- glm(formula = f_bal_identified_orem, data = lds_replication, family=binomial(link = "logit"))
summary(bal_identified_orem_mod)
# Salt Lake City
# specify formula
f_bal_identified_salt_lake_city <- formula(paste(var_bal_identified_loc[3], paste(c(var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
# logit model
bal_identified_salt_lake_city_mod <- glm(formula = f_bal_identified_salt_lake_city, data = lds_replication, family=binomial(link = "logit"))
summary(bal_identified_salt_lake_city_mod)
# Utah State
# specify formula
f_bal_identified_utah <- formula(paste(var_bal_identified_loc[4], paste(c(var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
# logit model
bal_identified_utah_mod <- glm(formula = f_bal_identified_utah, data = lds_replication, family=binomial(link = "logit"))
summary(bal_identified_utah_mod)
## Table A15: Demographics and Most Identified Places
# variable names
var_name_bal_identified_place_outcome <- c("Provo", "Orem", "Salt Lake City", "Utah State")
## NOTE: variable name is slightly different from the original draft ##
var_name_bal_identified_place <- var_name_demog
# regression table
stargazer(bal_identified_provo_mod,
          bal_identified_orem_mod, 
          bal_identified_salt_lake_city_mod, 
          bal_identified_utah_mod,
          omit = c("mission_year"),
          title = "Table A15: Demographics and Most Identified Places",
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          font.size = "small",
          covariate.labels = var_name_bal_identified_place, 
          dep.var.labels = var_name_bal_identified_place_outcome,
          add.lines = list(c("Mission Year FE?", "Yes", "Yes", "Yes")),
          label = "out_bal_identified_place")

## make Figure A3a: % Foreign-Born by State, 2010
# specify the data to create a choropleth map
acs_foreign_born_state_2010$value = acs_foreign_born_state_2010$est_foreign_born_2010  # set the desired 'value' column for choroplethr
acs_foreign_born_state_2010$region = tolower(acs_foreign_born_state_2010$state_name) # need lowercase state names
acs_foreign_born_state_2010_out <- acs_foreign_born_state_2010[c("value", "region")]
## Figure A3a: % Foreign-Born by State, 2010
state_choropleth(acs_foreign_born_state_2010_out,
                 title = "% Foreign-Born by State, 2010",
                 legend = NULL) +
  theme(plot.title = element_text(hjust = 0.5))

## make Figure A3b: % Latino by State, 2010
# specify the data to create a choropleth map
acs_latino_state_2010$value = acs_latino_state_2010$est_latino_2010  # set the desired 'value' column for choroplethr
acs_latino_state_2010$region = tolower(acs_latino_state_2010$state_name) # need lowercase state names
acs_latino_state_2010_out <- acs_latino_state_2010[c("value", "region")]
## Figure A3b: % Latino by State, 2010
state_choropleth(acs_latino_state_2010_out,
                 title = "% Latino by State, 2010",
                 legend = NULL) +
  theme(plot.title = element_text(hjust = 0.5))

## make Figure A4: State-Level Estimates of Immigration Attitudes, 2010-2011
# specify the data to create a choropleth map
cces_state_immi_est$value = cces_state_immi_est$est_state_immi_attitude  # set the desired 'value' column for choroplethr
cces_state_immi_est$region = tolower(cces_state_immi_est$state_name) # need lowercase state names
cces_state_immi_est_out <- cces_state_immi_est[c("value", "region")]
## Figure A4: State-Level Estimates of Immigration Attitudes, 2010-2011
state_choropleth(cces_state_immi_est_out,
                 title = "State-Level Estimates of Immigration Attitudes, 2010-2011",
                 legend = NULL) +
  theme(plot.title = element_text(hjust = 0.5))

## make Figure A5: State-Level Estimates of Immigration Attitudes, 2010-2011
# readjust the presentation order so that the most anti-immigrant state comes at the bottom
cces_state_immi_est <- cces_state_immi_est[order(cces_state_immi_est$est_state_immi_attitude), ]
cces_state_immi_est$presentation_order <- rep(NA)
cces_state_immi_est$presentation_order <- c(51:1)
## Figure A5: State-Level Estimates of Immigration Attitudes, 2010-2011
plot_cces_immi_attitude_state <- ggplot() + 
  geom_point(data = cces_state_immi_est, aes(x = est_state_immi_attitude, y = presentation_order), color = 'black') +
  geom_point(data = cces_state_immi_est[cces_state_immi_est$state_abb == "UT", ], aes(x = est_state_immi_attitude, y = presentation_order), size = 5) + # highlight Utah
  scale_y_continuous(breaks=1:51, labels = as.vector(rev(cces_state_immi_est$state_abb))) +
  ylab("State") +
  xlab("Anti-Immigrant Sentiment") +
  xlim(c(0, 1)) +
  ggtitle("State-Level Estimates of Immigration Attitudes, 2010-2011", subtitle = NULL) +
  theme_light() +
  theme(text = element_text(size = 16), 
        plot.title = element_text(hjust = 0.5))
print(plot_cces_immi_attitude_state)

## make Figure A6: 2012-2014 CCES Panel Study
# republican respondents
gop_cces_panel <- subset(cces_panel_2012_2014, cces_panel_2012_2014$gop == 1)
# LDS respondents
lds_cces_panel <- subset(cces_panel_2012_2014, cces_panel_2012_2014$lds == 1)
# mean anti-immigrant attitudes
# republican
gop_mean_immi_attitude_2012 <- mean(gop_cces_panel$immi_attitude_2012, na.rm = TRUE)
gop_mean_immi_attitude_2014 <- mean(gop_cces_panel$immi_attitude_2014, na.rm = TRUE)
# LDS
lds_mean_immi_attitude_2012 <- mean(lds_cces_panel$immi_attitude_2012, na.rm = TRUE)
lds_mean_immi_attitude_2014 <- mean(lds_cces_panel$immi_attitude_2014, na.rm = TRUE)
# standard errors
# republican
gop_se_mean_immi_attitude_2012 <- (sd(gop_cces_panel$immi_attitude_2012, na.rm = TRUE)) / sqrt(dim(gop_cces_panel[complete.cases(gop_cces_panel$immi_attitude_2012), ])[1])
gop_se_mean_immi_attitude_2014 <- (sd(gop_cces_panel$immi_attitude_2014, na.rm = TRUE)) / sqrt(dim(gop_cces_panel[complete.cases(gop_cces_panel$immi_attitude_2014), ])[1])
# LDS
lds_se_mean_immi_attitude_2012 <- (sd(lds_cces_panel$immi_attitude_2012, na.rm = TRUE)) / sqrt(dim(lds_cces_panel[complete.cases(lds_cces_panel$immi_attitude_2012), ])[1])
lds_se_mean_immi_attitude_2014 <- (sd(lds_cces_panel$immi_attitude_2014, na.rm = TRUE)) / sqrt(dim(lds_cces_panel[complete.cases(lds_cces_panel$immi_attitude_2014), ])[1])
# confidence intervals
# republican
gop_ci_mean_immi_attitude_2012 <- gop_mean_immi_attitude_2012 + (c(-1, 1) * qnorm(0.975) * gop_se_mean_immi_attitude_2012)
gop_ci_mean_immi_attitude_2014 <- gop_mean_immi_attitude_2014 + (c(-1, 1) * qnorm(0.975) * gop_se_mean_immi_attitude_2014)
# LDS
lds_ci_mean_immi_attitude_2012 <- lds_mean_immi_attitude_2012 + (c(-1, 1) * qnorm(0.975) * lds_se_mean_immi_attitude_2012)
lds_ci_mean_immi_attitude_2014 <- lds_mean_immi_attitude_2014 + (c(-1, 1) * qnorm(0.975) * lds_se_mean_immi_attitude_2014)
# combine results
# republican
vec_est_gop_2012 <- c(gop_mean_immi_attitude_2012, gop_se_mean_immi_attitude_2012, gop_ci_mean_immi_attitude_2012)
vec_est_gop_2014 <- c(gop_mean_immi_attitude_2014, gop_se_mean_immi_attitude_2014, gop_ci_mean_immi_attitude_2014)
# LDS
vec_est_lds_2012 <- c(lds_mean_immi_attitude_2012, lds_se_mean_immi_attitude_2012, lds_ci_mean_immi_attitude_2012)
vec_est_lds_2014 <- c(lds_mean_immi_attitude_2014, lds_se_mean_immi_attitude_2014, lds_ci_mean_immi_attitude_2014)
# matrix
mat_est_cces_panel <- rbind(vec_est_gop_2012,
                            vec_est_gop_2014,
                            vec_est_lds_2012,
                            vec_est_lds_2014)
# data frame
dat_est_cces_panel <- as.data.frame(mat_est_cces_panel)
row.names(dat_est_cces_panel) <- c()
# column names
colnames(dat_est_cces_panel) <- c("mean", "se", "ci_lower_95", "ci_upper_95")
# year
dat_est_cces_panel$year <- rep(NA)
dat_est_cces_panel$year <- as.vector(rep(c(2012, 2014), 2))
# respondent attributes
dat_est_cces_panel$`Respondent Attributes` <- rep(NA)
dat_est_cces_panel$`Respondent Attributes` <- as.vector(c(rep("Republicans", 2), rep("LDS", 2)))
## Figure A6: 2012-2014 CCES Panel Study
plot_cces_panel_2012_2014 <- ggplot() +
  geom_pointrange(data = dat_est_cces_panel, 
                  aes(x = year, 
                      y = mean,
                      ymin = ci_lower_95, 
                      ymax = ci_upper_95,
                      shape = `Respondent Attributes`),
                  size = 0.8) +
  geom_line(data = dat_est_cces_panel, 
            aes(x = year, 
                y = mean,
                lty = `Respondent Attributes`),
            size = 0.6) +
  ylab("Immigration Attitudes \n(Higher Values = More Hostile toward Immigrants)") +
  xlab("Year") +
  ylim(c(0, 1)) +
  geom_hline(aes(yintercept=0.5), linetype="dashed") +
  ggtitle("2012-2014 CCES Panel Study", subtitle = NULL) +
  scale_colour_brewer(palette = "Set1", direction = -1) +
  theme_light(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5))
print(plot_cces_panel_2012_2014)
## make Table A16a: Principal Component Analysis of Dependent Variable Items: Pre-Mission
# specify pre-mission attitudes for pca
var_immi_attitudes_pre_pca <- c("immi_english_pre",
                                "immi_tuition_pre",
                                "immi_welfare_pre",
                                "immi_employer_pre",
                                "immi_legal_pre",
                                "immi_violence_pre",
                                "immi_lazy_pre",
                                "immi_honest_pre",
                                "immi_victim_pre",
                                "immi_rights_pre")
# specify post-mission attitudes for pca
var_immi_attitudes_post_pca <- c("immi_english_post",
                                 "immi_tuition_post",
                                 "immi_welfare_post",
                                 "immi_employer_post",
                                 "immi_legal_post",
                                 "immi_violence_post",
                                 "immi_lazy_post",
                                 "immi_honest_post",
                                 "immi_victim_post",
                                 "immi_rights_post")
# select variables of interest
lds_replication_pre_pca <- lds_replication[var_immi_attitudes_pre_pca]
# drop missing data
lds_replication_pre_pca <- lds_replication_pre_pca[complete.cases(lds_replication_pre_pca), ]
# we first use principal function here because it allows us to prespecify the number of principal components
# pre-mission immigration attitudes
pca_out_pre <- principal(lds_replication_pre_pca, nfactors = 2, rotate="varimax")
# combine loadings
mat_pca_out_pre <- cbind(pca_out_pre$loadings[, 1], pca_out_pre$loadings[, 2])
dat_pca_out_pre <- as.data.frame(mat_pca_out_pre)
# column names
colnames(dat_pca_out_pre) <- c("First Principal Component", "Second Principal Component")
# row names
row.names(dat_pca_out_pre) <- c("English", "Tuition", "Welfare", "Employer", "Legal Status", "Violence", "Lazy", "Honest", "Victim", "Rights")
## Table A16a: Principal Component Analysis of Dependent Variable Items: Pre-Mission
# latex outputs
tab_pca_out_pre <- xtable(dat_pca_out_pre, digits = 2, caption = "Principal Component Analysis of Dependent Variable Items: Pre-Mission")
print(tab_pca_out_pre)

## make Table A16b: Principal Component Analysis of Dependent Variable Items: Post-Mission
# select variables of interest
lds_replication_post_pca <- lds_replication[var_immi_attitudes_post_pca]
# drop missing data
lds_replication_post_pca <- lds_replication_post_pca[complete.cases(lds_replication_post_pca), ]
# we first use principal function here because it allows us to prespecify the number of principal components
# post-mission immigration attitudes
pca_out_post <- principal(lds_replication_post_pca, nfactors = 2, rotate="varimax")
# combine loadings
mat_pca_out_post <- cbind(pca_out_post$loadings[, 1], pca_out_post$loadings[, 2])
dat_pca_out_post <- as.data.frame(mat_pca_out_post)
# column names
colnames(dat_pca_out_post) <- c("First Principal Component", "Second Principal Component")
# row names
row.names(dat_pca_out_post) <- c("English", "Tuition", "Welfare", "Employer", "Legal Status", "Violence", "Lazy", "Honest", "Victim", "Rights")
## Table A16b: Principal Component Analysis of Dependent Variable Items: Post-Mission
# latex outputs
tab_pca_out_post <- xtable(dat_pca_out_post, digits = 2, caption = "Principal Component Analysis of Dependent Variable Items: Post-Mission")
print(tab_pca_out_post)

## make Figure A7a: Scree Plot of Pre-Mission Immigration Items
# use prcomp function to produce scree plot: pre-mission attitudes
pca_out_pre_scree <- prcomp(lds_replication_pre_pca, scale = TRUE)
# scree plot
## Figure A7a: Scree Plot of Pre-Mission Immigration Items
fviz_eig(pca_out_pre_scree) +
  labs(title = "Scree Plot: Pre-Mission Items") +
  theme_light(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5))

## make Figure A7b: Scree Plot of Post-Mission Immigration Items
# use prcomp function to produce scree plot: post-mission attitudes
pca_out_post_scree <- prcomp(lds_replication_post_pca, scale = TRUE)
# scree plot
# figure a7b: scree plot of post-mission immigration items
fviz_eig(pca_out_post_scree) +
  labs(title = "Scree Plot: Post-Mission Items") +
  theme_light(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5))

# cronbach's alpha on page 18
# pre-mission policy items
var_immi_attitude_pre_policy <- c("immi_english_pre",
                                   "immi_tuition_pre",
                                   "immi_welfare_pre",
                                   "immi_legal_pre",
                                   "immi_employer_pre",
                                   "immi_rights_pre")
# select variables
dat_immi_attitude_pre_policy <- lds_replication[var_immi_attitude_pre_policy]
# cronbach's alpha
cronbach_alpha_immi_attitude_pre_policy <- alpha(dat_immi_attitude_pre_policy)
cronbach_alpha_immi_attitude_pre_policy <- cronbach_alpha_immi_attitude_pre_policy$total
cronbach_alpha_immi_attitude_pre_policy <- cronbach_alpha_immi_attitude_pre_policy$raw_alpha
print(cronbach_alpha_immi_attitude_pre_policy)
# 0.76 for pre-mission policy items

# post-mission policy items
var_immi_attitude_post_policy <- c("immi_english_post",
                                  "immi_tuition_post",
                                  "immi_welfare_post",
                                  "immi_legal_post",
                                  "immi_employer_post",
                                  "immi_rights_post")
# select variables
dat_immi_attitude_post_policy <- lds_replication[var_immi_attitude_post_policy]
# cronbach's alpha
cronbach_alpha_immi_attitude_post_policy <- alpha(dat_immi_attitude_post_policy)
cronbach_alpha_immi_attitude_post_policy <- cronbach_alpha_immi_attitude_post_policy$total
cronbach_alpha_immi_attitude_post_policy <- cronbach_alpha_immi_attitude_post_policy$raw_alpha
print(cronbach_alpha_immi_attitude_post_policy)
# 0.77 for post-mission policy items

# pre-mission stereotypes items
var_immi_attitude_pre_stereotype <- c("immi_violence_pre",
                                      "immi_lazy_pre",
                                      "immi_honest_pre")
# select variables
dat_immi_attitude_pre_stereotype <- lds_replication[var_immi_attitude_pre_stereotype]
# cronbach's alpha
cronbach_alpha_immi_attitude_pre_stereotype <- alpha(dat_immi_attitude_pre_stereotype)
cronbach_alpha_immi_attitude_pre_stereotype <- cronbach_alpha_immi_attitude_pre_stereotype$total
cronbach_alpha_immi_attitude_pre_stereotype <- cronbach_alpha_immi_attitude_pre_stereotype$raw_alpha
print(cronbach_alpha_immi_attitude_pre_stereotype)
# 0.59 for pre-mission stereotype items

# post-mission stereotypes items
var_immi_attitude_post_stereotype <- c("immi_violence_post",
                                       "immi_lazy_post",
                                       "immi_honest_post")
# select variables
dat_immi_attitude_post_stereotype <- lds_replication[var_immi_attitude_post_stereotype]
# cronbach's alpha
cronbach_alpha_immi_attitude_post_stereotype <- alpha(dat_immi_attitude_post_stereotype)
cronbach_alpha_immi_attitude_post_stereotype <- cronbach_alpha_immi_attitude_post_stereotype$total
cronbach_alpha_immi_attitude_post_stereotype <- cronbach_alpha_immi_attitude_post_stereotype$raw_alpha
print(cronbach_alpha_immi_attitude_post_stereotype)
# 0.60 for post-mission stereotype items

## make Figure A8: Distribution of Policy Items Index and Distribution of Stereotype Items Index
# descriptive statistics for plot
# statistics below also correspond to the statistics reported on page 19
# policy items
# mean
mean_immi_policy_pre <- mean(lds_replication$immi_policy_pre, na.rm = TRUE)
print(mean_immi_policy_pre)
mean_immi_policy_post <- mean(lds_replication$immi_policy_post, na.rm = TRUE)
print(mean_immi_policy_post)
# standard deviation
sd_immi_policy_pre <- sd(lds_replication$immi_policy_pre, na.rm = TRUE)
print(sd_immi_policy_pre)
sd_immi_policy_post <- sd(lds_replication$immi_policy_post, na.rm = TRUE)
print(sd_immi_policy_post)
# distribution of policy items: pre-mission
policy_pre_dist <- lds_replication$immi_policy_pre
policy_pre_dist <- as.data.frame(policy_pre_dist)
# pre/post indicator
policy_pre_dist$pre_post <- rep(NA)
policy_pre_dist$pre_post <- "Pre-Mission"
# add mean
policy_pre_dist$immi_policy_pre_mean <- rep(NA)
policy_pre_dist$immi_policy_pre_mean <- rep(mean_immi_policy_pre, dim(policy_pre_dist)[1])
# column names
colnames(policy_pre_dist) <- c("immi_policy_pre", "pre_post", "mean")
# distribution of policy items: post-mission
policy_post_dist <- lds_replication$immi_policy_post
policy_post_dist <- as.data.frame(policy_post_dist)
# pre/post indicator
policy_post_dist$pre_post <- rep(NA)
policy_post_dist$pre_post <- "Post-Mission"
# add mean
policy_post_dist$immi_policy_post_mean <- rep(NA)
policy_post_dist$immi_policy_post_mean <- rep(mean_immi_policy_post, dim(policy_post_dist)[1])
# column names
colnames(policy_post_dist) <- c("immi_policy_post", "pre_post", "mean")
# stereotype items
# mean
mean_immi_stereotype_pre <- mean(lds_replication$immi_stereotype_pre, na.rm = TRUE)
print(mean_immi_stereotype_pre)
mean_immi_stereotype_post <- mean(lds_replication$immi_stereotype_post, na.rm = TRUE)
print(mean_immi_stereotype_post)
# standard deviation
sd_immi_stereotype_pre <- sd(lds_replication$immi_stereotype_pre, na.rm = TRUE)
print(sd_immi_stereotype_pre)
sd_immi_stereotype_post <- sd(lds_replication$immi_stereotype_post, na.rm = TRUE)
print(sd_immi_stereotype_post)
# distribution of stereotype items: pre-mission
stereotype_pre_dist <- lds_replication$immi_stereotype_pre
stereotype_pre_dist <- as.data.frame(stereotype_pre_dist)
# pre/post indicator
stereotype_pre_dist$pre_post <- rep(NA)
stereotype_pre_dist$pre_post <- "Pre-Mission"
# add mean
stereotype_pre_dist$immi_stereotype_pre_mean <- rep(NA)
stereotype_pre_dist$immi_stereotype_pre_mean <- rep(mean_immi_stereotype_pre, dim(stereotype_pre_dist)[1])
# column names
colnames(stereotype_pre_dist) <- c("immi_stereotype_pre", "pre_post", "mean")
# distribution of stereotype items: post-mission
stereotype_post_dist <- lds_replication$immi_stereotype_post
stereotype_post_dist <- as.data.frame(stereotype_post_dist)
# pre/post indicator
stereotype_post_dist$pre_post <- rep(NA)
stereotype_post_dist$pre_post <- "Post-Mission"
# add mean
stereotype_post_dist$immi_stereotype_post_mean <- rep(NA)
stereotype_post_dist$immi_stereotype_post_mean <- rep(mean_immi_stereotype_post, dim(stereotype_post_dist)[1])
# column names
colnames(stereotype_post_dist) <- c("immi_stereotype_post", "pre_post", "mean")
## Figure A8: Distribution of Policy Items Index and Distribution of Stereotype Items Index
# pre-mission, policy items
plot_immi_policy_pre_dist <- ggplot(policy_pre_dist) +
  geom_histogram(aes(immi_policy_pre), color = "black", fill = "red", size = 1.5)  +
  geom_vline(xintercept = mean_immi_policy_pre, color = "black", linetype="dashed", size = 1.5) +
  labs(title = "Policy Items, Pre-Mission") +
  ylim(c(0, 700)) +
  xlim(c(0, 1)) +
  xlab(NULL) +
  ylab(NULL) +
  scale_x_continuous(breaks = seq(0, 1, by = 0.25)) +
  scale_x_continuous() +
  theme_light(base_size = 14) +
  theme(text = element_text(size = 16), 
        plot.title = element_text(hjust = 0.5))
# post-mission, policy items
plot_immi_policy_post_dist <- ggplot(policy_post_dist) +
  geom_histogram(aes(immi_policy_post), color = "black", fill = "blue", size = 1.5)  +
  geom_vline(xintercept = mean_immi_policy_post, color = "black", linetype="dashed", size = 1.5) +
  labs(title = "Policy Items, Post-Mission") +
  ylim(c(0, 700)) +
  xlim(c(0, 1)) +
  xlab(NULL) +
  ylab(NULL) +
  scale_x_continuous(breaks = seq(0, 1, by = 0.25)) +
  theme_light(base_size = 14) +
  theme(text = element_text(size = 16), 
        plot.title = element_text(hjust = 0.5))
# pre-mission, stereotype items
######################NOTE##########################
#### stereotype items have fewer unique values  ####
#### bar chart looks better and                 ####
#### substantively the same as histogram        ####
######################NOTE##########################
plot_immi_stereotypes_pre_dist <- ggplot(stereotype_pre_dist) +
  geom_bar(aes(immi_stereotype_pre), color = "black", fill = "red", size = 1.5)  +
  geom_vline(xintercept = mean_immi_stereotype_pre, color = "black", linetype="dashed", size = 1.5) +
  labs(title = "Stereotype Items, Pre-Mission") +
  xlab(NULL) +
  ylab(NULL) +
  ylim(c(0, 700)) +
  xlim(c(-0.1, 1)) +
  scale_x_continuous(breaks = seq(-0.1, 1, by = 0.25)) +
  theme_light(base_size = 14) +
  theme(text = element_text(size = 16), 
        plot.title = element_text(hjust = 0.5))

# post-mission, stereotype items
plot_immi_stereotypes_post_dist <- ggplot(stereotype_post_dist) +
  geom_bar(aes(immi_stereotype_post), color = "black", fill = "blue", size = 1.5)  +
  geom_vline(xintercept = mean_immi_stereotype_post, color = "black", linetype="dashed", size = 1.5) +
  labs(title = "Stereotype Items, Post-Mission") +
  xlab(NULL) +
  ylab(NULL) +
  ylim(c(0, 700)) +
  xlim(c(-0.1, 1)) +
  theme_light(base_size = 14) +
  theme(text = element_text(size = 16), 
        plot.title = element_text(hjust = 0.5))
# combine plots
bind_immi_attitude_dist <- ggarrange(plot_immi_policy_pre_dist, 
                                     plot_immi_policy_post_dist,
                                     plot_immi_stereotypes_pre_dist,
                                     plot_immi_stereotypes_post_dist,
                                     widths = c(0.3, 0.3),
                                     heights = c(0.5, 0.5), 
                                     ncol = 2, 
                                     nrow = 2)

## Figure A8: Distribution of Policy Items Index and Distribution of Stereotype Items Index
plot_immi_attitude_dist <- annotate_figure(bind_immi_attitude_dist,
                                         top = text_grob("Distributions of Dependent Variables", 
                                                         color = "Black", 
                                                         face = "bold", 
                                                         size = 20),
                                         bottom = text_grob("Immigration Attitudes", 
                                                            color = "Black", 
                                                            face = "bold", 
                                                            size = 18),
                                         left = text_grob("Frequency", color = "Black", face = "bold", size = 16, rot = 90))
print(plot_immi_attitude_dist)

## make Figure A9: Immigration and Anti???Immigrant Attitudes
## Figure A9: Immigration and Anti???Immigrant Attitudes
# use country_state_county_desc data frame
plot_desc_immi <- ggplot() + 
  geom_smooth(data = country_state_county_desc, 
              aes(x = foreign_born_2010, y = mean, color = Group, lty = Group), 
              size = 1.5, 
              method = "lm",
              formula = y ~ x, se = FALSE) +
  ylab("Anti-Immigrant Attitudes") +
  xlab("% Foreign Born Population") +
  ylim(c(0, 1)) +
  ggtitle("Immigration and Anti-Immigrant Attitudes", subtitle = NULL) +
  theme_light() +
  theme(text = element_text(size = 16), 
        plot.title = element_text(hjust = 0.5))
print(plot_desc_immi)

## make Table A17: Pre-Post Comparisons
# need to convert the data to long version for pre-post comparison
lds_replication_pre <- as.data.frame(lds_replication)
lds_replication_post <- as.data.frame(lds_replication)
# create a immigration attitude variable that can be combined
# policy attitudes
# pre-mission
lds_replication_pre$immi_policy <- rep(NA)
lds_replication_pre$immi_policy <- as.vector(as.numeric(lds_replication_pre$immi_policy_pre))
# post-mission
lds_replication_post$immi_policy <- rep(NA)
lds_replication_post$immi_policy <- as.vector(as.numeric(lds_replication_post$immi_policy_post))
# stereotype attitudes
# pre-mission
lds_replication_pre$immi_stereotype <- rep(NA)
lds_replication_pre$immi_stereotype <- as.vector(as.numeric(lds_replication_pre$immi_stereotype_pre))
# post-mission
lds_replication_post$immi_stereotype <- rep(NA)
lds_replication_post$immi_stereotype <- as.vector(as.numeric(lds_replication_post$immi_stereotype_post))
# add a variable indicating post survey 
lds_replication_pre$post <- rep(NA)
lds_replication_pre$post <- rep(0, dim(lds_replication_pre)[1])
lds_replication_post$post <- rep(NA)
lds_replication_post$post <- rep(1, dim(lds_replication_post)[1])
# create the long version
lds_replication_pre_post_long <- rbind(lds_replication_pre, lds_replication_post)
dim(lds_replication_pre_post_long)
# usa sample
lds_replication_pre_post_long_usa <- subset(lds_replication_pre_post_long, lds_replication_pre_post_long$usa == 1)
# abroad sample
lds_replication_pre_post_long_abroad <- subset(lds_replication_pre_post_long, lds_replication_pre_post_long$usa == 0)
# simple pre-post comparison: full sample, policy items
# specify formula
f_pre_post_policy_full <- formula(paste("immi_policy", paste(c("post", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_policy_full_mod <- lm(formula = f_pre_post_policy_full, data = lds_replication_pre_post_long)
summary(pre_post_policy_full_mod)
# simple pre-post comparison: full sample, stereotype items
# specify formula
f_pre_post_stereotype_full <- formula(paste("immi_stereotype", paste(c("post", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_stereotype_full_mod <- lm(formula = f_pre_post_stereotype_full, data = lds_replication_pre_post_long)
summary(pre_post_stereotype_full_mod)
# simple pre-post comparison: USA sample, policy items
# specify formula
f_pre_post_policy_usa <- formula(paste("immi_policy", paste(c("post", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_policy_usa_mod <- lm(formula = f_pre_post_policy_usa, data = lds_replication_pre_post_long_usa)
summary(pre_post_policy_usa_mod)
# simple pre-post comparison: USA sample, stereotype items
# specify formula
f_pre_post_stereotype_usa <- formula(paste("immi_stereotype", paste(c("post", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_stereotype_usa_mod <- lm(formula = f_pre_post_stereotype_usa, data = lds_replication_pre_post_long_usa)
summary(pre_post_stereotype_usa_mod)
# simple pre-post comparison: Abroad sample, policy items
# specify formula
f_pre_post_policy_abroad <- formula(paste("immi_policy", paste(c("post", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_policy_abroad_mod <- lm(formula = f_pre_post_policy_abroad, data = lds_replication_pre_post_long_abroad)
summary(pre_post_policy_abroad_mod)
# simple pre-post comparison: Abroad sample, stereotype items
# specify formula
f_pre_post_stereotype_abroad <- formula(paste("immi_stereotype", paste(c("post", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_stereotype_abroad_mod <- lm(formula = f_pre_post_stereotype_abroad, data = lds_replication_pre_post_long_abroad)
summary(pre_post_stereotype_abroad_mod)
## Table A17: Pre-Post Comparisons
# variable names
var_name_pre_post_outcome <- c("Policy/Full", "Stereotype/Full", "Policy/USA", "Stereotype/USA", "Policy/Abroad", "Stereotype/Abroad")
var_name_pre_post <- c("Post", var_name_demog)
# regression table
# long version inflates the number of observations. I manually added the n
stargazer(pre_post_policy_full_mod, 
          pre_post_stereotype_full_mod, 
          pre_post_policy_usa_mod, 
          pre_post_stereotype_usa_mod, 
          pre_post_policy_abroad_mod, 
          pre_post_stereotype_abroad_mod,
          omit = c("mission_year"),
          title = "Pre-Post Comparisons",
          digits = 3,
          keep.stat = c("n"), 
          omit.stat = c("n"), 
          style = "apsr", 
          font.size = "small",
          covariate.labels = var_name_pre_post, 
          dep.var.labels = var_name_pre_post_outcome,
          add.lines = list(c("Mission Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("N", "1402", "1400", "468", "467", "882", "881")), 
          label = "pre_post_comparisons")

## make Table A18: DID Estimates: Effects of Serving in USA
# policy items
# specify formula
# no covariate
f_usa_abroad_policy_1 <- formula(paste("diff_immi_policy", paste(c("usa"), collapse = " + "), sep = " ~ "))
usa_abroad_policy_mod_1 <- lm(formula = f_usa_abroad_policy_1, data = lds_replication)
summary(usa_abroad_policy_mod_1)
# with covariates
f_usa_abroad_policy_2 <- formula(paste("diff_immi_policy", paste(c("usa", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
usa_abroad_policy_mod_2 <- lm(formula = f_usa_abroad_policy_2, data = lds_replication)
summary(usa_abroad_policy_mod_2)
# stereotype items
# specify formula
# no covariate
f_usa_abroad_stereotype_1 <- formula(paste("diff_immi_stereotype", paste(c("usa"), collapse = " + "), sep = " ~ "))
usa_abroad_stereotype_mod_1 <- lm(formula = f_usa_abroad_stereotype_1, data = lds_replication)
summary(usa_abroad_stereotype_mod_1)
# with covariates
f_usa_abroad_stereotype_2 <- formula(paste("diff_immi_stereotype", paste(c("usa", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
usa_abroad_stereotype_mod_2 <- lm(formula = f_usa_abroad_stereotype_2, data = lds_replication)
summary(usa_abroad_stereotype_mod_2)

## Table A18: DID Estimates: Effects of Serving in USA
# variable names
var_usa_abroad_effect_outcome <- c("Policy", "Stereotype")
var_usa_abroad_effect <- c("Served in USA", var_name_demog)
# regression table
stargazer(usa_abroad_policy_mod_1,
          usa_abroad_policy_mod_2,
          usa_abroad_stereotype_mod_1,
          usa_abroad_stereotype_mod_2,
          title = "DID Estimates: Effects of Serving in USA",
          omit = c("mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_usa_abroad_effect, 
          dep.var.labels = var_usa_abroad_effect_outcome,
          font.size = "small",
          add.lines = list(c("Year FE?", "No", "Yes", "No", "Yes")),
          label = "usa_abroad_effect")

## make Figure A10: DID estimates: Abroad vs. US.
# extract coefficients to create plot
# policy items
# no covariate
# point estimate
coef_usa_abroad_policy_1 <- summary(usa_abroad_policy_mod_1)$coefficients[2, 1]
# standard error
se_usa_abroad_policy_1<- summary(usa_abroad_policy_mod_1)$coefficients[2, 2]
# combine
qoi_usa_abroad_policy_1 <- as.vector(c(coef_usa_abroad_policy_1, se_usa_abroad_policy_1))
# confidence intervals
ci_usa_abroad_policy_1  <- as.vector(confint(usa_abroad_policy_mod_1)[2, ])
# add confidence intervals
qoi_usa_abroad_policy_1 <- c(qoi_usa_abroad_policy_1, 
                             ci_usa_abroad_policy_1)
qoi_usa_abroad_policy_1 <- as.data.frame(qoi_usa_abroad_policy_1)
# readjust
qoi_usa_abroad_policy_1 <- t(qoi_usa_abroad_policy_1)
qoi_usa_abroad_policy_1 <- as.data.frame(qoi_usa_abroad_policy_1)
# column names
colnames(qoi_usa_abroad_policy_1) <- c("point_est", 
                                       "se", 
                                       "usa_abroad_ci_lower_95", 
                                       "usa_abroad_ci_upper_95")
# clean row names
row.names(qoi_usa_abroad_policy_1) <- c()
# policy items
# with covariates
# point estimate
coef_usa_abroad_policy_2 <- summary(usa_abroad_policy_mod_2)$coefficients[2, 1]
# standard error
se_usa_abroad_policy_2<- summary(usa_abroad_policy_mod_2)$coefficients[2, 2]
# combine
qoi_usa_abroad_policy_2 <- as.vector(c(coef_usa_abroad_policy_2, se_usa_abroad_policy_2))
# confidence intervals
ci_usa_abroad_policy_2  <- as.vector(confint(usa_abroad_policy_mod_2)[2, ])
# add confidence intervals
qoi_usa_abroad_policy_2 <- c(qoi_usa_abroad_policy_2, 
                             ci_usa_abroad_policy_2)
qoi_usa_abroad_policy_2 <- as.data.frame(qoi_usa_abroad_policy_2)
# readjust
qoi_usa_abroad_policy_2 <- t(qoi_usa_abroad_policy_2)
qoi_usa_abroad_policy_2 <- as.data.frame(qoi_usa_abroad_policy_2)
# column names
colnames(qoi_usa_abroad_policy_2) <- c("point_est", 
                                       "se", 
                                       "usa_abroad_ci_lower_95", 
                                       "usa_abroad_ci_upper_95")
# clean row names
row.names(qoi_usa_abroad_policy_2) <- c()
# stereotype items
# no covariate
# point estimate
coef_usa_abroad_stereotype_1 <- summary(usa_abroad_stereotype_mod_1)$coefficients[2, 1]
# standard error
se_usa_abroad_stereotype_1<- summary(usa_abroad_stereotype_mod_1)$coefficients[2, 2]
# combine
qoi_usa_abroad_stereotype_1 <- as.vector(c(coef_usa_abroad_stereotype_1, se_usa_abroad_stereotype_1))
# confidence intervals
ci_usa_abroad_stereotype_1  <- as.vector(confint(usa_abroad_stereotype_mod_1)[2, ])
# add confidence intervals
qoi_usa_abroad_stereotype_1 <- c(qoi_usa_abroad_stereotype_1, 
                             ci_usa_abroad_stereotype_1)
qoi_usa_abroad_stereotype_1 <- as.data.frame(qoi_usa_abroad_stereotype_1)
# readjust
qoi_usa_abroad_stereotype_1 <- t(qoi_usa_abroad_stereotype_1)
qoi_usa_abroad_stereotype_1 <- as.data.frame(qoi_usa_abroad_stereotype_1)
# column names
colnames(qoi_usa_abroad_stereotype_1) <- c("point_est", 
                                       "se", 
                                       "usa_abroad_ci_lower_95", 
                                       "usa_abroad_ci_upper_95")
# clean row names
row.names(qoi_usa_abroad_stereotype_1) <- c()
# stereotype items
# with covariates
# point estimate
coef_usa_abroad_stereotype_2 <- summary(usa_abroad_stereotype_mod_2)$coefficients[2, 1]
# standard error
se_usa_abroad_stereotype_2<- summary(usa_abroad_stereotype_mod_2)$coefficients[2, 2]
# combine
qoi_usa_abroad_stereotype_2 <- as.vector(c(coef_usa_abroad_stereotype_2, se_usa_abroad_stereotype_2))
# confidence intervals
ci_usa_abroad_stereotype_2  <- as.vector(confint(usa_abroad_stereotype_mod_2)[2, ])
# add confidence intervals
qoi_usa_abroad_stereotype_2 <- c(qoi_usa_abroad_stereotype_2, 
                             ci_usa_abroad_stereotype_2)
qoi_usa_abroad_stereotype_2 <- as.data.frame(qoi_usa_abroad_stereotype_2)
# readjust
qoi_usa_abroad_stereotype_2 <- t(qoi_usa_abroad_stereotype_2)
qoi_usa_abroad_stereotype_2 <- as.data.frame(qoi_usa_abroad_stereotype_2)
# column names
colnames(qoi_usa_abroad_stereotype_2) <- c("point_est", 
                                       "se", 
                                       "usa_abroad_ci_lower_95", 
                                       "usa_abroad_ci_upper_95")
# clean row names
row.names(qoi_usa_abroad_stereotype_2) <- c()
# combine results
qoi_usa_abroad <- rbind(qoi_usa_abroad_policy_1,
                        qoi_usa_abroad_policy_2,
                        qoi_usa_abroad_stereotype_1,
                        qoi_usa_abroad_stereotype_2)
qoi_usa_abroad <- as.data.frame(qoi_usa_abroad)
# add domain indicator
qoi_usa_abroad$domain <- rep(NA)
qoi_usa_abroad$domain <- c(rep("Policy", 2), rep("Stereotype", 2))
# covariate indicator
qoi_usa_abroad$covariates <- rep(NA)
qoi_usa_abroad$covariates <- c(rep(c("No Covariate", "With Covariates"), 2))
# subset by policy and stereotype
qoi_usa_abroad_policy <- subset(qoi_usa_abroad, qoi_usa_abroad$domain == "Policy")
qoi_usa_abroad_stereotype <- subset(qoi_usa_abroad, qoi_usa_abroad$domain == "Stereotype")
## Figure A10: DID estimates: Abroad vs. US.
# policy
# annotated text of point estimates
qoi_usa_abroad_policy$point_est_text <- rep(NA)
qoi_usa_abroad_policy$point_est_text <- as.character(round(qoi_usa_abroad_policy$point_est, 3))
# policy item plot
plot_qoi_usa_abroad_policy <- ggplot() +
  geom_pointrange(data = qoi_usa_abroad_policy, 
                  aes(x = covariates, 
                      y = point_est, 
                      ymin = usa_abroad_ci_lower_95, 
                      ymax = usa_abroad_ci_upper_95,
                      shape = covariates),
                  position = position_dodge(width = 0.4),
                  size = 0.8) +
  annotate("text",
           x = c(2.1, 0.9), # adjust position
           y = c(-0.066, -0.056), # adjust position
           label = c(rev(qoi_usa_abroad_policy$point_est_text))) +
  ylab(NULL) +
  xlab(NULL) +
  geom_hline(yintercept=0, linetype="dashed") +
  ylim(c(-0.15, 0.15)) +
  ggtitle("Policy", subtitle = NULL) +
  theme_light(base_size = 14) +
  coord_flip() +
  theme(text = element_text(size = 16), 
        plot.title = element_text(hjust = 0.5))
# stereotype
# annotated text of point estimates
qoi_usa_abroad_stereotype$point_est_text <- rep(NA)
qoi_usa_abroad_stereotype$point_est_text <- as.character(round(qoi_usa_abroad_stereotype$point_est, 3))
# stereotype item plot
plot_qoi_usa_abroad_stereotype <- ggplot() +
  geom_pointrange(data = qoi_usa_abroad_stereotype, 
                  aes(x = covariates, 
                      y = point_est, 
                      ymin = usa_abroad_ci_lower_95, 
                      ymax = usa_abroad_ci_upper_95,
                      shape = covariates),
                  position = position_dodge(width = 0.4),
                  size = 0.8) +
  annotate("text",
           x = c(2.1, 0.9), # adjust position
           y = c(-0.060, -0.050), # adjust position
           label = c(rev(qoi_usa_abroad_stereotype$point_est_text))) +
  ylab(NULL) +
  xlab(NULL) +
  geom_hline(yintercept=0, linetype="dashed") +
  ylim(c(-0.15, 0.15)) +
  ggtitle("Stereotype", subtitle = NULL) +
  theme_light(base_size = 14) +
  coord_flip() +
  theme(axis.text.y = element_blank(),
        text = element_text(size = 16), 
        plot.title = element_text(hjust = 0.5))

# combine plots
bind_plot_usa_abroad_mod <- ggarrange(plot_qoi_usa_abroad_policy, 
                                      plot_qoi_usa_abroad_stereotype,
                                      widths = c(1, 1),
                                      ncol = 2, 
                                      nrow = 1,
                                      common.legend = TRUE, 
                                      align = "hv",
                                      legend = "right")

# Figure A10: DID Estimates: Abroad vs. US.
out_plot_usa_abroad_mod <- annotate_figure(bind_plot_usa_abroad_mod,
                                           top = text_grob("Effects of Serving in USA", 
                                                           color = "Black", 
                                                           face = "bold", 
                                                           size = 20),
                                                    bottom = text_grob("Change in Immigration Attitudes", 
                                                                       color = "Black", 
                                                                       face = "bold", 
                                                                       size = 18),
                                                    left = text_grob("Covariates", color = "Black", face = "bold", size = 16, rot = 90))
print(out_plot_usa_abroad_mod)

## make Figure A11: Effects of serving in the USA: Individual Survey Items
############################################################
# individual survey items
var_individual_item <- c("diff_immi_english",
                         "diff_immi_tuition",
                         "diff_immi_welfare",
                         "diff_immi_employer",
                         "diff_immi_legal",
                         "diff_immi_violence",
                         "diff_immi_lazy",
                         "diff_immi_honest",
                         "diff_immi_rights")
# select variables above
y_bind_individual_item <- lds_replication[var_individual_item]
# to matrix
y_bind_individual_item <- as.matrix(y_bind_individual_item)
# individual item
individual_item_usa_abroad_mod <- lm(y_bind_individual_item ~ 
                                     usa +
                                     female +
                                     age +
                                     white +
                                     family_income +
                                     prior_health_issue +
                                     prior_other_lang +
                                     pid_republican +
                                     ideo_conservative +
                                     as.factor(mission_year),
                                    data = lds_replication)
sum_individual_item_usa_abroad_mod <- summary(individual_item_usa_abroad_mod)
# confidence intervals
ci_individual_item <- confint(individual_item_usa_abroad_mod)
ci_individual_item <- as.data.frame(ci_individual_item)
# subset by row names
ci_individual_item$row_name <- rep(NA)
ci_individual_item$row_name <- as.vector(as.character(row.names(ci_individual_item)))
# confidence intervals on usa
ci_individual_item_usa <- subset(ci_individual_item, ci_individual_item$row_name %in%
                                   paste(var_individual_item, "usa", sep = ":"))
# extract point estimates
point_est_individual_item_usa <- coefficients(individual_item_usa_abroad_mod)[2, ]
point_est_individual_item_usa <- as.vector(as.numeric(point_est_individual_item_usa))
# combine point estimates and confidence intervals
mat_est_individual_item_usa <- cbind(point_est_individual_item_usa,
                                     ci_individual_item_usa$`2.5 %`,
                                     ci_individual_item_usa$`97.5 %`)
dat_est_individual_item_usa <- as.data.frame(mat_est_individual_item_usa)
colnames(dat_est_individual_item_usa) <- c("point_est", "ci_lower_95", "ci_upper_95")
# add names of survey items for plot
var_name_individual_item <- factor(c("English",
                                     "Tuition",
                                     "Welfare",
                                     "Employer",
                                     "Legal Status",
                                     "Violence",
                                     "Lazy",
                                     "Honest",
                                     "Rights"), 
                                   levels = rev(c("English",
                                                   "Tuition",
                                                   "Welfare",
                                                   "Employer",
                                                   "Legal Status",
                                                   "Violence",
                                                   "Lazy",
                                                   "Honest",
                                                   "Rights")))
# add to the data frame
dat_est_individual_item_usa$item_name <- rep(NA)
dat_est_individual_item_usa$item_name <- var_name_individual_item
# text version of point estimate for annotation
dat_est_individual_item_usa$point_est_text <- rep(NA)
dat_est_individual_item_usa$point_est_text <- as.character(round(dat_est_individual_item_usa$point_est, 3))
# Figure A11: Effects of serving in the USA: Individual Survey Items
plot_usa_abroad_individual_item <- ggplot() +
  geom_pointrange(data = dat_est_individual_item_usa, 
                  aes(x = item_name, 
                      y = point_est, 
                      ymin = ci_lower_95, 
                      ymax = ci_upper_95),
                  position = position_dodge(width = 0.4),
                  size = 0.8) +
  annotate("text",
           x = c(9.28, 8.28, 7.28, 6.28, 5.28, 4.28, 3.28, 2.28, 1.28), # specify the position of annotation
           y = dat_est_individual_item_usa$point_est,
           label = dat_est_individual_item_usa$point_est_text) +
  ylab("Change in Immigration Attitudes") +
  xlab("Estimates") +
  geom_hline(yintercept = 0, linetype="dashed") +
  ylim(c(-0.2, 0.2)) +
  ggtitle("Effects of serving in the USA: Individual Survey Items") +
  scale_colour_brewer(palette = "Set1", direction=-1) +
  theme_light(base_size = 14) +
  coord_flip() +
  theme(text = element_text(size = 16), 
        plot.title = element_text(hjust = 0.5))
print(plot_usa_abroad_individual_item)

## make Table A19: Effects of Larger Immigrant Communities
# this is essentially the full version of Table 4
## Table A19: Effects of Larger Immigrant Communities
# variable names
var_name_intl_migrant_outcome_full <- c("Policy", "Stereotype")
var_name_intl_migrant_full <- c("% International Migrant", var_name_demog)
# regression table
stargazer(intl_migrant_policy_mod_1,
          intl_migrant_policy_mod_2,
          intl_migrant_stereotype_mod_1,
          intl_migrant_stereotype_mod_2,
          title = "Effects of Larger Immigrant Communities (% International Migrant)",
          omit = c("mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_intl_migrant_full, 
          dep.var.labels = var_name_intl_migrant_outcome_full,
          font.size = "small",
          add.lines = list(c("Demographics", "No", "Yes", "No", "Yes"), 
                           c("Year FE?", "No", "Yes", "No", "Yes")),
          label = "diff_imm_intl_migrant_full")

## make Table A20: Effects of Serving in Latin American Countries
# policy items
# specify formula
# no covariate
f_latin_america_policy_1 <- formula(paste("diff_immi_policy", paste(c("latin_america"), collapse = " + "), sep = " ~ "))
latin_america_policy_mod_1 <- lm(formula = f_latin_america_policy_1, data = lds_replication)
summary(latin_america_policy_mod_1)
# with covariates
f_latin_america_policy_2 <- formula(paste("diff_immi_policy", paste(c("latin_america", "non_latin_america", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
latin_america_policy_mod_2 <- lm(formula = f_latin_america_policy_2, data = lds_replication)
summary(latin_america_policy_mod_2)
# stereotype items
# specify formula
# no covariate
f_latin_america_stereotype_1 <- formula(paste("diff_immi_stereotype", paste(c("latin_america"), collapse = " + "), sep = " ~ "))
latin_america_stereotype_mod_1 <- lm(formula = f_latin_america_stereotype_1, data = lds_replication)
summary(latin_america_stereotype_mod_1)
# with covariates
f_latin_america_stereotype_2 <- formula(paste("diff_immi_stereotype", paste(c("latin_america", "non_latin_america", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
latin_america_stereotype_mod_2 <- lm(formula = f_latin_america_stereotype_2, data = lds_replication)
summary(latin_america_stereotype_mod_2)

## Table A20: Effects of Serving in Latin American Countries
# variable names
var_name_latin_america_outcome <- c("Policy", "Stereotype")
var_name_latin_america <- c("Latin America", "Non-Latin America", var_name_demog)
# regression table
stargazer(latin_america_policy_mod_1, 
          latin_america_policy_mod_2, 
          latin_america_stereotype_mod_1, 
          latin_america_stereotype_mod_2,
          omit = c("mission_year"),
          title = "Effects of Serving in Latin American Countries",
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_latin_america, 
          dep.var.labels = var_name_latin_america_outcome,
          font.size = "small",
          add.lines = list(c("Year FE?", "No", "Yes", "No", "Yes")),
          label = "latin_america")

## statistics on page 24
# non-english assignment
table(lds_replication_usa$non_eng_mission, exclude = NULL)
# 136 out of 575
print(table(lds_replication_usa$non_eng_mission, exclude = NULL)[2] / nrow(lds_replication_usa)) * 100
# roughly 24%
# among those assigned to speak a foreign language
lds_replication_usa_non_eng_mission <- subset(lds_replication_usa, lds_replication_usa$non_eng_mission == 1)
table(lds_replication_usa_non_eng_mission$esp_mission, exclude = NULL)
print(table(lds_replication_usa_non_eng_mission$esp_mission, exclude = NULL)[2] / nrow(lds_replication_usa_non_eng_mission) * 100)
# 87.5% were assigned to speak spanish
# among those assigned to speak spanish
lds_replication_usa_esp_mission <- subset(lds_replication_usa, lds_replication_usa$esp_mission == 1)
table(lds_replication_usa_esp_mission$speak_esp, exclude = NULL)
# only one respondent ended up not speaking spanish

## make Table A21: Serving in US and Non-English Language Assignment
# specify formula
f_non_eng_usa <- formula(paste("non_eng_mission", paste(c("usa", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
# logit model
non_eng_usa_mod <- glm(formula = f_non_eng_usa, data = lds_replication, family=binomial(link = "logit"))
summary(non_eng_usa_mod)
## Table A21: Serving in US and Non-English Language Assignment
# variable names
var_name_non_eng_usa_outcome <- c("Non-English Assignment")
var_name_non_eng_usa <- c("Served in USA", var_name_demog)
# regression table
stargazer(non_eng_usa_mod,
          title = "Serving Abroad and Speaking Foreign Languages",
          omit = c("mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_non_eng_usa, 
          dep.var.labels = var_name_non_eng_usa_outcome,
          font.size = "small",
          add.lines = list(c("Year FE?", "Yes")),
          label = "non_eng_usa")

## make Table A22: Effects of being Assigned to Speak a Foreign Language
# already have the results for abroad sample
# focus on those served abroad
# policy items
# specify formula
# with covariates
f_non_eng_mission_abroad_policy <- formula(paste(var_diff_immi_outcome[1], paste(c("non_eng_mission", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_abroad_policy_mod <- lm(formula = f_non_eng_mission_abroad_policy, data = lds_replication_abroad)
summary(non_eng_mission_abroad_policy_mod)
# stereotype items
# specify formula
# with covariates
f_non_eng_mission_abroad_stereotype <- formula(paste(var_diff_immi_outcome[2], paste(c("non_eng_mission", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_abroad_stereotype_mod <- lm(formula = f_non_eng_mission_abroad_stereotype, data = lds_replication_abroad)
summary(non_eng_mission_abroad_stereotype_mod)
## Table A22: Effects of being Assigned to Speak a Foreign Language
# variable names
var_name_non_eng_mission_outcome <- c("Policy/USA", "Stereotype/USA", "Policy/Abroad", "Stereotype/Abroad")
var_name_non_eng_mission <- c("Assigned to Speak a Foreign Language", var_name_demog)
# regression table
stargazer(non_eng_mission_usa_policy_mod,
          non_eng_mission_usa_stereotype_mod,
          non_eng_mission_abroad_policy_mod,
          non_eng_mission_abroad_stereotype_mod,
          title = "Effects of being Assigned to Speak a Foreign Language",
          omit = c("mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_non_eng_mission, 
          dep.var.labels = var_name_non_eng_mission_outcome,
          font.size = "small",
          add.lines = list(c("Year FE?", "Yes", "Yes", "Yes", "Yes")),
          label = "non_eng_lang")

## make Table A23: Effects of Serving in USA on Spanish Speaker
# respondents assigned to speak spanish
lds_replication_esp_mission <- subset(lds_replication, lds_replication$esp_mission == 1)
# policy items
# specify formula
# with covariates
f_esp_mission_policy <- formula(paste(var_diff_immi_outcome[1], paste(c("usa", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
esp_mission_policy_mod <- lm(formula = f_esp_mission_policy, data = lds_replication_esp_mission)
summary(esp_mission_policy_mod)
# stereotype items
# specify formula
# with covariates
f_esp_mission_stereotype <- formula(paste(var_diff_immi_outcome[2], paste(c("usa", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
esp_mission_stereotype_mod <- lm(formula = f_esp_mission_stereotype, data = lds_replication_esp_mission)
summary(esp_mission_stereotype_mod)
# respondents who atucally spoke spanish while serving the mission
lds_replication_speak_esp <- subset(lds_replication, lds_replication$speak_esp == 1)
# policy items
# specify formula
# with covariates
f_speak_esp_policy <- formula(paste(var_diff_immi_outcome[1], paste(c("usa", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
speak_esp_policy_mod <- lm(formula = f_speak_esp_policy, data = lds_replication_speak_esp)
summary(speak_esp_policy_mod)
# stereotype items
# specify formula
# with covariates
f_speak_esp_stereotype <- formula(paste(var_diff_immi_outcome[2], paste(c("usa", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
speak_esp_stereotype_mod <- lm(formula = f_speak_esp_stereotype, data = lds_replication_speak_esp)
summary(speak_esp_stereotype_mod)
## Table A23: Effects of Serving in USA on Spanish Speaker
# variable names
var_name_usa_lang_esp_outcome <- c("Policy/Assign", "Stereotype/Assign", "Policy/Speak", "Stereotype/Speak")
var_name_usa_lang_esp <- c("Served in USA", var_name_demog)
# regression table
stargazer(esp_mission_policy_mod,
          esp_mission_stereotype_mod,
          speak_esp_policy_mod,
          speak_esp_stereotype_mod,
          title = "Effects of Serving in USA on Spanish Speaker",
          omit = c("mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_usa_lang_esp, 
          dep.var.labels = var_name_usa_lang_esp_outcome,
          font.size = "small",
          add.lines = list(c("Year FE?", "Yes", "Yes", "Yes", "Yes")),
          label = "usa_lang_esp")

## make Figure A12: Placebo Test: Serving in USA and Abortion Attitudes
# specify formula
# with covariates
f_usa_abortion <- formula(paste("diff_abortion", paste(c("usa", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
usa_abortion_mod <- lm(formula = f_usa_abortion, data = lds_replication)
summary(usa_abortion_mod)
# extract coefficient
coef_usa_abortion <- summary(usa_abortion_mod)$coefficients[2, 1]
se_usa_abortion <- summary(usa_abortion_mod)$coefficients[2, 2]
# combine point estimates and standard errors
qoi_usa_abortion <- as.vector(c(coef_usa_abortion, se_usa_abortion))
# confidence intervals
ci_usa_abotion <- confint(usa_abortion_mod)[2, ]
# combine point estimates, standard errors, and confidence intervals
qoi_usa_abortion <- c(qoi_usa_abortion, 
                      ci_usa_abotion)
qoi_usa_abortion <- as.data.frame(qoi_usa_abortion)
# readjust
qoi_usa_abortion <- t(qoi_usa_abortion)
qoi_usa_abortion <- as.data.frame(qoi_usa_abortion)
# column names
colnames(qoi_usa_abortion) <- c("point_est", 
                                "se", 
                                "ci_usa_abortion_lower_95", 
                                "ci_usa_abortion_upper_95")
# variable name
qoi_usa_abortion$var_name_abortion <- rep(NA)
qoi_usa_abortion$var_name_abortion <- "Abortion"
# clean up row name
row.names(qoi_usa_abortion) <- c()
## Figure A12: Placebo Test: Serving in USA and Abortion Attitudes
plot_lds_usa_abortion <- ggplot() +
  geom_pointrange(data = qoi_usa_abortion, 
                  aes(x = var_name_abortion, 
                      y = point_est, 
                      ymin = ci_usa_abortion_lower_95, 
                      ymax = ci_usa_abortion_upper_95),
                  position = position_dodge(width = 0.4),
                  size = 0.8) +
  ylab("Coefficient") +
  xlab("Issue") +
  geom_hline(yintercept=0, linetype="dashed") +
  ylim(c(-0.2, 0.2)) +
  ggtitle("Serving in USA and Abortion Attitudes", subtitle = NULL) +
  theme_light(base_size = 14) +
  coord_flip() +
  theme(plot.title = element_text(hjust = 0.5))
print(plot_lds_usa_abortion)

## make Figure A13: Placebo Test: Speaking Spanish and Abortion Attitudes. Only US participants included
# specify formula
# with covariates
f_usa_esp_mission_abortion <- formula(paste("diff_abortion", paste(c("esp_mission", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
usa_esp_mission_abortion_mod <- lm(formula = f_usa_esp_mission_abortion, data = lds_replication_usa)
summary(usa_esp_mission_abortion_mod)
# extract coefficient
coef_usa_esp_mission_abortion <- summary(usa_esp_mission_abortion_mod)$coefficients[2, 1]
se_usa_esp_mission_abortion <- summary(usa_esp_mission_abortion_mod)$coefficients[2, 2]
# combine point estimates and standard errors
qoi_usa_esp_mission_abortion <- as.vector(c(coef_usa_esp_mission_abortion, se_usa_esp_mission_abortion))
# confidence intervals
ci_usa_esp_mission_abotion <- confint(usa_esp_mission_abortion_mod)[2, ]
# combine point estimates, standard errors, and confidence intervals
qoi_usa_esp_mission_abortion <- c(qoi_usa_esp_mission_abortion, 
                      ci_usa_esp_mission_abotion)
qoi_usa_esp_mission_abortion <- as.data.frame(qoi_usa_esp_mission_abortion)
# readjust
qoi_usa_esp_mission_abortion <- t(qoi_usa_esp_mission_abortion)
qoi_usa_esp_mission_abortion <- as.data.frame(qoi_usa_esp_mission_abortion)
# column names
colnames(qoi_usa_esp_mission_abortion) <- c("point_est", 
                                            "se", 
                                            "ci_usa_esp_mission_abortion_lower_95", 
                                            "ci_usa_esp_mission_abortion_upper_95")
# variable name
qoi_usa_esp_mission_abortion$var_name_abortion <- rep(NA)
qoi_usa_esp_mission_abortion$var_name_abortion <- "Abortion"
# clean up row name
row.names(qoi_usa_esp_mission_abortion) <- c()

## Figure A13: Placebo Test: Speaking Spanish and Abortion Attitudes. Only US participants included
plot_lds_usa_esp_mission_abortion <- ggplot() +
  geom_pointrange(data = qoi_usa_esp_mission_abortion, 
                  aes(x = var_name_abortion, 
                      y = point_est, 
                      ymin = ci_usa_esp_mission_abortion_lower_95, 
                      ymax = ci_usa_esp_mission_abortion_upper_95),
                  position = position_dodge(width = 0.4),
                  size = 0.8) +
  ylab("Coefficient") +
  xlab("Issue") +
  geom_hline(yintercept=0, linetype="dashed") +
  ylim(c(-0.2, 0.2)) +
  ggtitle("Speaking Spanish and Abortion Attitudes", subtitle = NULL) +
  theme_light(base_size = 14) +
  coord_flip() +
  theme(plot.title = element_text(hjust = 0.5))
print(plot_lds_usa_esp_mission_abortion)

# make Table A24a: Pre-Post Comparisons: Male Missionaries
# gender subset
# male missionaries
lds_replication_male <- subset(lds_replication, lds_replication$female == 0)
# usa and abroad
lds_replication_usa_male <- subset(lds_replication_male, lds_replication_male$usa == 1)
lds_replication_abroad_male <- subset(lds_replication_male, lds_replication_male$usa == 0)
# female missionaries
lds_replication_female <- subset(lds_replication, lds_replication$female == 1)
# usa and abroad
lds_replication_usa_female <- subset(lds_replication_female, lds_replication_female$usa == 1)
lds_replication_abroad_female <- subset(lds_replication_female, lds_replication_female$usa == 0)
# male subset
# need to convert the data to long version for pre-post comparison
lds_replication_pre_male <- as.data.frame(lds_replication_male)
lds_replication_post_male <- as.data.frame(lds_replication_male)
# create a immigration attitude variable that can be combined
# policy attitudes
# pre-mission
lds_replication_pre_male$immi_policy <- rep(NA)
lds_replication_pre_male$immi_policy <- as.vector(as.numeric(lds_replication_pre_male$immi_policy_pre))
# post-mission
lds_replication_post_male$immi_policy <- rep(NA)
lds_replication_post_male$immi_policy <- as.vector(as.numeric(lds_replication_post_male$immi_policy_post))
# stereotype attitudes
# pre-mission
lds_replication_pre_male$immi_stereotype <- rep(NA)
lds_replication_pre_male$immi_stereotype <- as.vector(as.numeric(lds_replication_pre_male$immi_stereotype_pre))
# post-mission
lds_replication_post_male$immi_stereotype <- rep(NA)
lds_replication_post_male$immi_stereotype <- as.vector(as.numeric(lds_replication_post_male$immi_stereotype_post))
# add a variable indicating post survey 
lds_replication_pre_male$post <- rep(NA)
lds_replication_pre_male$post <- rep(0, dim(lds_replication_pre_male)[1])
lds_replication_post_male$post <- rep(NA)
lds_replication_post_male$post <- rep(1, dim(lds_replication_post_male)[1])
# create the long version
lds_replication_pre_post_long_male <- rbind(lds_replication_pre_male, lds_replication_post_male)
dim(lds_replication_pre_post_long_male)
# usa sample
lds_replication_pre_post_long_usa_male <- subset(lds_replication_pre_post_long_male, lds_replication_pre_post_long_male$usa == 1)
# abroad sample
lds_replication_pre_post_long_abroad_male <- subset(lds_replication_pre_post_long_male, lds_replication_pre_post_long_male$usa == 0)
# simple pre-post comparison: full sample, policy items
# specify formula: drop gender as a control
f_pre_post_policy_full_male <- formula(paste("immi_policy", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_policy_full_mod_male <- lm(formula = f_pre_post_policy_full_male, data = lds_replication_pre_post_long_male)
summary(pre_post_policy_full_mod_male)
# simple pre-post comparison: full sample, stereotype items
# specify formula
f_pre_post_stereotype_full_male <- formula(paste("immi_stereotype", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_stereotype_full_mod_male <- lm(formula = f_pre_post_stereotype_full_male, data = lds_replication_pre_post_long_male)
summary(pre_post_stereotype_full_mod_male)
# simple pre-post comparison: USA sample, policy items
# specify formula
f_pre_post_policy_usa_male <- formula(paste("immi_policy", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_policy_usa_mod_male <- lm(formula = f_pre_post_policy_usa_male, data = lds_replication_pre_post_long_usa_male)
summary(pre_post_policy_usa_mod_male)
# simple pre-post comparison: USA sample, stereotype items
# specify formula
f_pre_post_stereotype_usa_male <- formula(paste("immi_stereotype", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_stereotype_usa_mod_male <- lm(formula = f_pre_post_stereotype_usa_male, data = lds_replication_pre_post_long_usa_male)
summary(pre_post_stereotype_usa_mod_male)
# simple pre-post comparison: Abroad sample, policy items
# specify formula
f_pre_post_policy_abroad_male <- formula(paste("immi_policy", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_policy_abroad_mod_male <- lm(formula = f_pre_post_policy_abroad_male, data = lds_replication_pre_post_long_abroad_male)
summary(pre_post_policy_abroad_mod_male)
# simple pre-post comparison: Abroad sample, stereotype items
# specify formula
f_pre_post_stereotype_abroad_male <- formula(paste("immi_stereotype", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_stereotype_abroad_mod_male <- lm(formula = f_pre_post_stereotype_abroad_male, data = lds_replication_pre_post_long_abroad_male)
summary(pre_post_stereotype_abroad_mod_male)
## Table A24a: Pre-Post Comparisons: Male Missionaries
# variable names
var_name_pre_post_outcome_male <- c("Policy/Full", "Stereotype/Full", "Policy/USA", "Stereotype/USA", "Policy/Abroad", "Stereotype/Abroad")
var_name_pre_post_male <- c("Post", var_name_demog[! var_name_demog %in% c("Female")])
# regression table
# long version inflates the number of observations. I manually added the n
stargazer(pre_post_policy_full_mod_male, 
          pre_post_stereotype_full_mod_male, 
          pre_post_policy_usa_mod_male, 
          pre_post_stereotype_usa_mod_male, 
          pre_post_policy_abroad_mod_male, 
          pre_post_stereotype_abroad_mod_male,
          omit = c("mission_year"),
          title = "Pre-Post Comparisons: Male Missionaries",
          digits = 3,
          keep.stat = c("n"), 
          omit.stat = c("n"), 
          style = "apsr", 
          font.size = "small",
          covariate.labels = var_name_pre_post_male, 
          dep.var.labels = var_name_pre_post_outcome_male,
          add.lines = list(c("Mission Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("N", "417", "416", "89", "89", "307", "306")), 
          label = "pre_post_comparisons_male")

## make Table A24a: Pre-Post Comparisons: Female Missionaries
# female subset
# need to convert the data to long version for pre-post comparison
lds_replication_pre_female <- as.data.frame(lds_replication_female)
lds_replication_post_female <- as.data.frame(lds_replication_female)
# create a immigration attitude variable that can be combined
# policy attitudes
# pre-mission
lds_replication_pre_female$immi_policy <- rep(NA)
lds_replication_pre_female$immi_policy <- as.vector(as.numeric(lds_replication_pre_female$immi_policy_pre))
# post-mission
lds_replication_post_female$immi_policy <- rep(NA)
lds_replication_post_female$immi_policy <- as.vector(as.numeric(lds_replication_post_female$immi_policy_post))
# stereotype attitudes
# pre-mission
lds_replication_pre_female$immi_stereotype <- rep(NA)
lds_replication_pre_female$immi_stereotype <- as.vector(as.numeric(lds_replication_pre_female$immi_stereotype_pre))
# post-mission
lds_replication_post_female$immi_stereotype <- rep(NA)
lds_replication_post_female$immi_stereotype <- as.vector(as.numeric(lds_replication_post_female$immi_stereotype_post))
# add a variable indicating post survey 
lds_replication_pre_female$post <- rep(NA)
lds_replication_pre_female$post <- rep(0, dim(lds_replication_pre_female)[1])
lds_replication_post_female$post <- rep(NA)
lds_replication_post_female$post <- rep(1, dim(lds_replication_post_female)[1])
# create the long version
lds_replication_pre_post_long_female <- rbind(lds_replication_pre_female, lds_replication_post_female)
dim(lds_replication_pre_post_long_female)
# usa sample
lds_replication_pre_post_long_usa_female <- subset(lds_replication_pre_post_long_female, lds_replication_pre_post_long_female$usa == 1)
# abroad sample
lds_replication_pre_post_long_abroad_female <- subset(lds_replication_pre_post_long_female, lds_replication_pre_post_long_female$usa == 0)
# simple pre-post comparison: full sample, policy items
# specify formula: drop gender as a control
f_pre_post_policy_full_female <- formula(paste("immi_policy", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_policy_full_mod_female <- lm(formula = f_pre_post_policy_full_female, data = lds_replication_pre_post_long_female)
summary(pre_post_policy_full_mod_female)
# simple pre-post comparison: full sample, stereotype items
# specify formula
f_pre_post_stereotype_full_female <- formula(paste("immi_stereotype", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_stereotype_full_mod_female <- lm(formula = f_pre_post_stereotype_full_female, data = lds_replication_pre_post_long_female)
summary(pre_post_stereotype_full_mod_female)
# simple pre-post comparison: USA sample, policy items
# specify formula
f_pre_post_policy_usa_female <- formula(paste("immi_policy", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_policy_usa_mod_female <- lm(formula = f_pre_post_policy_usa_female, data = lds_replication_pre_post_long_usa_female)
summary(pre_post_policy_usa_mod_female)
# simple pre-post comparison: USA sample, stereotype items
# specify formula
f_pre_post_stereotype_usa_female <- formula(paste("immi_stereotype", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_stereotype_usa_mod_female <- lm(formula = f_pre_post_stereotype_usa_female, data = lds_replication_pre_post_long_usa_female)
summary(pre_post_stereotype_usa_mod_female)
# simple pre-post comparison: Abroad sample, policy items
# specify formula
f_pre_post_policy_abroad_female <- formula(paste("immi_policy", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_policy_abroad_mod_female <- lm(formula = f_pre_post_policy_abroad_female, data = lds_replication_pre_post_long_abroad_female)
summary(pre_post_policy_abroad_mod_female)
# simple pre-post comparison: Abroad sample, stereotype items
# specify formula
f_pre_post_stereotype_abroad_female <- formula(paste("immi_stereotype", paste(c("post", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
pre_post_stereotype_abroad_mod_female <- lm(formula = f_pre_post_stereotype_abroad_female, data = lds_replication_pre_post_long_abroad_female)
summary(pre_post_stereotype_abroad_mod_female)
## Table A24a: Pre-Post Comparisons: Female Missionaries
# variable names
var_name_pre_post_outcome_female <- c("Policy/Full", "Stereotype/Full", "Policy/USA", "Stereotype/USA", "Policy/Abroad", "Stereotype/Abroad")
var_name_pre_post_female <- c("Post", var_name_demog[! var_name_demog %in% c("Female")])
# regression table
# long version inflates the number of observations. I manually added the n
stargazer(pre_post_policy_full_mod_female, 
          pre_post_stereotype_full_mod_female, 
          pre_post_policy_usa_mod_female, 
          pre_post_stereotype_usa_mod_female, 
          pre_post_policy_abroad_mod_female, 
          pre_post_stereotype_abroad_mod_female,
          omit = c("mission_year"),
          title = "Pre-Post Comparisons: Female Missionaries",
          digits = 3,
          keep.stat = c("n"), 
          omit.stat = c("n"), 
          style = "apsr", 
          font.size = "small",
          covariate.labels = var_name_pre_post_female, 
          dep.var.labels = var_name_pre_post_outcome_female,
          add.lines = list(c("Mission Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("N", "985", "984", "379", "378", "575", "574")), 
          label = "pre_post_comparisons_female")

## make Table A25: Effects of Serving in US: Split by Gender
# male missionaries
# policy items
# specify formula
# with covariates
f_usa_abroad_policy_male <- formula(paste("diff_immi_policy", paste(c("usa", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
usa_abroad_policy_mod_male <- lm(formula = f_usa_abroad_policy_male, data = lds_replication_male)
summary(usa_abroad_policy_mod_male)
# stereotype items
# specify formula
# with covariates
f_usa_abroad_stereotype_male <- formula(paste("diff_immi_stereotype", paste(c("usa", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
usa_abroad_stereotype_mod_male <- lm(formula = f_usa_abroad_stereotype_male, data = lds_replication_male)
summary(usa_abroad_stereotype_mod_male)
# female missionaries
# policy items
# specify formula
# with covariates
f_usa_abroad_policy_female <- formula(paste("diff_immi_policy", paste(c("usa", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
usa_abroad_policy_mod_female <- lm(formula = f_usa_abroad_policy_female, data = lds_replication_female)
summary(usa_abroad_policy_mod_female)
# stereotype items
# specify formula
# with covariates
f_usa_abroad_stereotype_female <- formula(paste("diff_immi_stereotype", paste(c("usa", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
usa_abroad_stereotype_mod_female <- lm(formula = f_usa_abroad_stereotype_female, data = lds_replication_female)
summary(usa_abroad_stereotype_mod_female)

## Table A25: Effects of Serving in US: Split by Gender
# variable names
var_name_usa_abroad_effect_outcome_gender <- c("Policy/Male", "Stereotype/Male", "Policy/Female", "Stereotype/Female")
var_name_usa_abroad_effect_gender <- c("Served in USA", var_name_demog[! var_name_demog %in% c("Female")])
# regression table
stargazer(usa_abroad_policy_mod_male,
          usa_abroad_stereotype_mod_male,
          usa_abroad_policy_mod_female,
          usa_abroad_stereotype_mod_female,
          title = "Effects of Serving in US: Split by Gender",
          omit = c("mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_usa_abroad_effect_gender, 
          dep.var.labels = var_name_usa_abroad_effect_outcome_gender,
          font.size = "small",
          add.lines = list(c("Year FE?", "Yes", "Yes", "Yes", "Yes")),
          label = "usa_abroad_effect_gender")

## make Table A26: Effects of Larger Immigrant Communities: Full Sample Split by Gender
# male subset
# policy items
# with covariates
f_intl_migrant_policy_male <- formula(paste(var_diff_immi_outcome[1], paste(c("loc_country_intl_migrant", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
intl_migrant_policy_mod_male <- lm(formula = f_intl_migrant_policy_male, data = lds_replication_male)
summary(intl_migrant_policy_mod_male)
# stereotype items
# with covariates
f_intl_migrant_stereotype_male <- formula(paste(var_diff_immi_outcome[2], paste(c("loc_country_intl_migrant", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
intl_migrant_stereotype_mod_male <- lm(formula = f_intl_migrant_stereotype_male, data = lds_replication_male)
summary(intl_migrant_stereotype_mod_male)
# female subset
# policy items
# with covariates
f_intl_migrant_policy_female <- formula(paste(var_diff_immi_outcome[1], paste(c("loc_country_intl_migrant", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
intl_migrant_policy_mod_female <- lm(formula = f_intl_migrant_policy_female, data = lds_replication_female)
summary(intl_migrant_policy_mod_female)
# stereotype items
# with covariates
f_intl_migrant_stereotype_female <- formula(paste(var_diff_immi_outcome[2], paste(c("loc_country_intl_migrant", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
intl_migrant_stereotype_mod_female <- lm(formula = f_intl_migrant_stereotype_female, data = lds_replication_female)
summary(intl_migrant_stereotype_mod_female)
## Table A26: Effects of Larger Immigrant Communities (International Migrant): Full Sample Split by Gender
# variable names
var_name_intl_migrant_outcome_gender <- c("Policy/Male", "Stereotype/Male", "Policy/Female", "Stereotype/Female")
var_name_intl_migrant_gender <- c("% International Migrant", var_name_demog[! var_name_demog %in% c("Female")])
# regression table
stargazer(intl_migrant_policy_mod_male,
          intl_migrant_stereotype_mod_male,
          intl_migrant_policy_mod_female,
          intl_migrant_stereotype_mod_female,
          title = "Effects of Larger Immigrant Communities (International Migrant): Full Sample Split by Gender",
          omit = c("mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_intl_migrant_gender, 
          dep.var.labels = var_name_intl_migrant_outcome_gender,
          font.size = "small",
          add.lines = list(c("Demographics", "Yes", "Yes", "Yes", "Yes"), 
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes")),
          label = "diff_imm_intl_migrant_gender")

## make Table A27a: Effects of Larger Latino/Immigrant Communities: Male Missionaries in US
# male subset
# specify formula
# latino, policy items
f_diff_county_latino_policy_male <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_loc_county_latino", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_latino_policy_mod_male <- lm(formula = f_diff_county_latino_policy_male, data = lds_replication_usa_male)
summary(diff_county_latino_policy_mod_male)
# foreign-born, policy items
f_diff_county_foreign_born_policy_male <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_loc_county_foreign_born", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_foreign_born_policy_mod_male <- lm(formula = f_diff_county_foreign_born_policy_male, data = lds_replication_usa_male)
summary(diff_county_foreign_born_policy_mod_male)
# foreign-born from latin american countrie, policy items
f_diff_county_fbla_policy_male <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_loc_county_fbla", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_fbla_policy_mod_male <- lm(formula = f_diff_county_fbla_policy_male, data = lds_replication_usa_male)
summary(diff_county_fbla_policy_mod_male)
# specify formula
# latino, stereotype items
f_diff_county_latino_stereotype_male <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_loc_county_latino", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_latino_stereotype_mod_male <- lm(formula = f_diff_county_latino_stereotype_male, data = lds_replication_usa_male)
summary(diff_county_latino_stereotype_mod_male)
# foreign-born, stereotype items
f_diff_county_foreign_born_stereotype_male <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_loc_county_foreign_born", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_foreign_born_stereotype_mod_male <- lm(formula = f_diff_county_foreign_born_stereotype_male, data = lds_replication_usa_male)
summary(diff_county_foreign_born_stereotype_mod_male)
# foreign-born from latin american countrie, stereotype items
f_diff_county_fbla_stereotype_male <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_loc_county_fbla", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_fbla_stereotype_mod_male <- lm(formula = f_diff_county_fbla_stereotype_male, data = lds_replication_usa_male)
summary(diff_county_fbla_stereotype_mod_male)
## Table A27a: Effects of Larger Latino/Immigrant Communities: Male Missionaries in US
# variable names
var_name_county_diff <- c("Diff. in % Latino", "Diff. in % Foreign Born", "Diff in % Foreign Born from Latin America")
var_name_county_diff_outcome_male <- rep(c("Policy", "Stereotype"), 3)
# regression table
stargazer(diff_county_latino_policy_mod_male, 
          diff_county_latino_stereotype_mod_male, 
          diff_county_foreign_born_policy_mod_male,
          diff_county_foreign_born_stereotype_mod_male,
          diff_county_fbla_policy_mod_male,
          diff_county_fbla_stereotype_mod_male,
          omit = c(var_demog[! var_demog %in% c("female")],
                   "mission_year"),
          title = "Effects of Larger Latino/Immigrant Communities: Male Missionaries in US",
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          font.size = "small",
          covariate.labels = var_name_county_diff, 
          dep.var.labels = var_name_county_diff_outcome_male,
          add.lines = list(c("Demographics?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), c("Mission Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
          label = "out_county_diff_mod_male")

## make Table A27b: Effects of Larger Latino/Immigrant Communities: Female Missionaries in US
# female subset
# specify formula
# latino, policy items
f_diff_county_latino_policy_female <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_loc_county_latino", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_latino_policy_mod_female <- lm(formula = f_diff_county_latino_policy_female, data = lds_replication_usa_female)
summary(diff_county_latino_policy_mod_female)
# foreign-born, policy items
f_diff_county_foreign_born_policy_female <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_loc_county_foreign_born", var_demog[! var_demog %in% c("fefemale")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_foreign_born_policy_mod_female <- lm(formula = f_diff_county_foreign_born_policy_female, data = lds_replication_usa_female)
summary(diff_county_foreign_born_policy_mod_female)
# foreign-born from latin american countrie, policy items
f_diff_county_fbla_policy_female <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_loc_county_fbla", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_fbla_policy_mod_female <- lm(formula = f_diff_county_fbla_policy_female, data = lds_replication_usa_female)
summary(diff_county_fbla_policy_mod_female)
# specify formula
# latino, stereotype items
f_diff_county_latino_stereotype_female <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_loc_county_latino", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_latino_stereotype_mod_female <- lm(formula = f_diff_county_latino_stereotype_female, data = lds_replication_usa_female)
summary(diff_county_latino_stereotype_mod_female)
# foreign-born, stereotype items
f_diff_county_foreign_born_stereotype_female <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_loc_county_foreign_born", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_foreign_born_stereotype_mod_female <- lm(formula = f_diff_county_foreign_born_stereotype_female, data = lds_replication_usa_female)
summary(diff_county_foreign_born_stereotype_mod_female)
# foreign-born from latin american countrie, stereotype items
f_diff_county_fbla_stereotype_female <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_loc_county_fbla", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
diff_county_fbla_stereotype_mod_female <- lm(formula = f_diff_county_fbla_stereotype_female, data = lds_replication_usa_female)
summary(diff_county_fbla_stereotype_mod_female)
## Table A27B: Effects of Larger Latino/Immigrant Communities: Female Missionaries in US
# variable names
var_name_county_diff <- c("Diff. in % Latino", "Diff. in % Foreign Born", "Diff in % Foreign Born from Latin America")
var_name_county_diff_outcome_female <- rep(c("Policy", "Stereotype"), 3)
# regression table
stargazer(diff_county_latino_policy_mod_female, 
          diff_county_latino_stereotype_mod_female, 
          diff_county_foreign_born_policy_mod_female,
          diff_county_foreign_born_stereotype_mod_female,
          diff_county_fbla_policy_mod_female,
          diff_county_fbla_stereotype_mod_female,
          omit = c(var_demog[! var_demog %in% c("Female")],
                   "mission_year"),
          title = "Effects of Larger Latino/Immigrant Communities: Female Missionaries in US",
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          font.size = "small",
          covariate.labels = var_name_county_diff, 
          dep.var.labels = var_name_county_diff_outcome_female,
          add.lines = list(c("Demographics?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), c("Mission Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
          label = "out_county_diff_mod_female")

## make Table A28a: Effects of being Assigned to Speak a Foreign Language: Male Missionaries
# run the regression model
# male missionaries in usa
# policy items
# with covariates
f_non_eng_mission_usa_policy_male <- formula(paste(var_diff_immi_outcome[1], paste(c("non_eng_mission", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_usa_policy_mod_male <- lm(formula = f_non_eng_mission_usa_policy_male, data = lds_replication_usa_male)
summary(non_eng_mission_usa_policy_mod_male)
# stereotype items
# with covariates
f_non_eng_mission_usa_stereotype_male <- formula(paste(var_diff_immi_outcome[2], paste(c("non_eng_mission", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_usa_stereotype_mod_male <- lm(formula = f_non_eng_mission_usa_stereotype_male, data = lds_replication_usa_male)
summary(non_eng_mission_usa_stereotype_mod_male)
# male missionaries abroad
# policy items
# with covariates
f_non_eng_mission_abroad_policy_male <- formula(paste(var_diff_immi_outcome[1], paste(c("non_eng_mission", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_abroad_policy_mod_male <- lm(formula = f_non_eng_mission_abroad_policy_male, data = lds_replication_abroad_male)
summary(non_eng_mission_abroad_policy_mod_male)
# stereotype items
# with covariates
f_non_eng_mission_abroad_stereotype_male <- formula(paste(var_diff_immi_outcome[2], paste(c("non_eng_mission", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_abroad_stereotype_mod_male <- lm(formula = f_non_eng_mission_abroad_stereotype_male, data = lds_replication_abroad_male)
summary(non_eng_mission_abroad_stereotype_mod_male)
## Table A28a: Effects of being Assigned to Speak a Foreign Language: Male Missionaries
# variable names
var_name_non_eng_mission_outcome_male <- c("Policy/USA", "Stereotype/USA", "Policy/Abroad", "Stereotype/Abroad")
var_name_non_eng_mission_male_male <- c("Assigned to Speak a Foreign Language", var_name_demog[! var_name_demog %in% c("Female")])
# regression table
stargazer(non_eng_mission_usa_policy_mod_male,
          non_eng_mission_usa_stereotype_mod_male,
          non_eng_mission_abroad_policy_mod_male,
          non_eng_mission_abroad_stereotype_mod_male,
          title = "Effects of being Assigned to Speak a Foreign Language: Male Missionaries",
          omit = c("mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_non_eng_mission_male_male, 
          dep.var.labels = var_name_non_eng_mission_outcome_male,
          font.size = "small",
          add.lines = list(c("Year FE?", "Yes", "Yes", "Yes", "Yes")),
          label = "non_eng_lang_male")

## make Table A28b: Effects of being Assigned to Speak a Foreign Language: Female Missionaries
# run the regression model
# female missionaries in usa
# policy items
# with covariates
f_non_eng_mission_usa_policy_female <- formula(paste(var_diff_immi_outcome[1], paste(c("non_eng_mission", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_usa_policy_mod_female <- lm(formula = f_non_eng_mission_usa_policy_female, data = lds_replication_usa_female)
summary(non_eng_mission_usa_policy_mod_female)
# stereotype items
# with covariates
f_non_eng_mission_usa_stereotype_female <- formula(paste(var_diff_immi_outcome[2], paste(c("non_eng_mission", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_usa_stereotype_mod_female <- lm(formula = f_non_eng_mission_usa_stereotype_female, data = lds_replication_usa_female)
summary(non_eng_mission_usa_stereotype_mod_female)
# female missionaries abroad
# policy items
# with covariates
f_non_eng_mission_abroad_policy_female <- formula(paste(var_diff_immi_outcome[1], paste(c("non_eng_mission", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_abroad_policy_mod_female <- lm(formula = f_non_eng_mission_abroad_policy_female, data = lds_replication_abroad_female)
summary(non_eng_mission_abroad_policy_mod_female)
# stereotype items
# with covariates
f_non_eng_mission_abroad_stereotype_female <- formula(paste(var_diff_immi_outcome[2], paste(c("non_eng_mission", var_demog[! var_demog %in% c("female")], "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
non_eng_mission_abroad_stereotype_mod_female <- lm(formula = f_non_eng_mission_abroad_stereotype_female, data = lds_replication_abroad_female)
summary(non_eng_mission_abroad_stereotype_mod_female)
## Table A28b: Effects of being Assigned to Speak a Foreign Language: Female Missionaries
# variable names
var_name_non_eng_mission_outcome_female <- c("Policy/USA", "Stereotype/USA", "Policy/Abroad", "Stereotype/Abroad")
var_name_non_eng_mission_female_female <- c("Assigned to Speak a Foreign Language", var_name_demog[! var_name_demog %in% c("Female")])
# regression table
stargazer(non_eng_mission_usa_policy_mod_female,
          non_eng_mission_usa_stereotype_mod_female,
          non_eng_mission_abroad_policy_mod_female,
          non_eng_mission_abroad_stereotype_mod_female,
          title = "Effects of being Assigned to Speak a Foreign Language: Female Missionaries",
          omit = c("mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_non_eng_mission_female_female, 
          dep.var.labels = var_name_non_eng_mission_outcome_female,
          font.size = "small",
          add.lines = list(c("Year FE?", "Yes", "Yes", "Yes", "Yes")),
          label = "non_eng_lang_female")

## make Table A29: Persistence of Effects: US Sample
# specify formula
# policy items
# foreign-born from Latin American countries
f_presist_fbla_policy <- formula(paste(var_diff_immi_outcome[1], paste(c("diff_loc_county_fbla * diff_month_complete", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
presist_fbla_policy_mod <- lm(formula = f_presist_fbla_policy, data = lds_replication_usa)
summary(presist_fbla_policy_mod)
# stereotype items
# foreign-born from Latin American countries
f_presist_fbla_stereotype <- formula(paste(var_diff_immi_outcome[2], paste(c("diff_loc_county_fbla * diff_month_complete", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
presist_fbla_stereotype_mod <- lm(formula = f_presist_fbla_stereotype, data = lds_replication_usa)
summary(presist_fbla_stereotype_mod)
# policy items
# assigned to speak a foreign language
f_presist_non_eng_mission_policy <- formula(paste(var_diff_immi_outcome[1], paste(c("non_eng_mission * diff_month_complete", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
presist_non_eng_mission_policy_mod <- lm(formula = f_presist_non_eng_mission_policy, data = lds_replication_usa)
summary(presist_non_eng_mission_policy_mod)
# stereotype items
# assigned to speak a foreign language
f_presist_non_eng_mission_stereotype <- formula(paste(var_diff_immi_outcome[2], paste(c("non_eng_mission * diff_month_complete", var_demog, "as.factor(mission_year)"), collapse = " + "), sep = " ~ "))
presist_non_eng_mission_stereotype_mod <- lm(formula = f_presist_non_eng_mission_stereotype, data = lds_replication_usa)
summary(presist_non_eng_mission_stereotype_mod)
## Table A29: Persistence of Effects: US Sample
# variable names
var_name_persist <- c("Diff. in % Foreign Born from Latin America", "Assigned to Speak a Foreign Language", "Diff. in % Months", "Diff. in % Foreign Born from Latin America * Diff. in % Months", "Assigned to Speak a Foreign Language * Diff. in % Months")
var_name_persist_outcome <- c("Policy", "Stereotype", "Policy", "Stereotype")
stargazer(presist_fbla_policy_mod,
          presist_fbla_stereotype_mod,
          presist_non_eng_mission_policy_mod,
          presist_non_eng_mission_stereotype_mod,
          title = "Persistence of Effects: US Sample",
          omit = c(var_demog, "mission_year"),
          digits = 3,
          keep.stat = c("n"), 
          style = "apsr", 
          covariate.labels = var_name_persist, 
          dep.var.labels = var_name_persist_outcome,
          font.size = "small",
          add.lines = list(c("Demographics", "Yes", "Yes", "Yes", "Yes"), 
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes")),
          label = "persist_effect")
